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Abstract We discuss the dynamics of zonal (or unidirectional) jets for barotropic flows forced by Gaussian stochastic fields 
with white in time correlation functions. This problem contains the stochastic dynamics of 2D Navier-Stokes equation as a 
special case. We consider the limit of weak forces and dissipation, when there is a time scale separation between the inertial time 
scale (fast) and the spin-up or spin-down time (large) needed to reach an average energy balance. In this limit, we show that an 
adiabatic reduction (or stochastic averaging) of the dynamics can be performed. We then obtain a kinetic equation that describes 
the slow evolution of zonal jets over a very long time scale, where the effect of non-zonal turbulence has been integrated out. 
\ The main theoretical difficulty, achieved in this work, is to analyze the stationary distribution of a Lyapunov equation that 
describes quasi-Gaussian fluctuations around each zonal jet, in the inertial limit. This is necessary to prove that there is no 
ultraviolet divergence at leading order in such a way that the asymptotic expansion is self-consistent. We obtain at leading order 
a Fokker-Planck equation, associated to a stochastic kinetic equation, that describes the slow jet dynamics. Its deterministic 
part is related to well known phenomenological theories (for instance Stochastic Structural Stability Theory) and to quasi-linear 
approximations, whereas the stochastic part allows to go beyond the computation of the most probable zonal jet. We argue that 
the effect of the stochastic part may be of huge importance when, as for instance in the proximity of phase transitions, more than 
one attractor of the dynamics is present. 
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•i-H ' 1 Introduction 



Turbulence in planetary atmospheres leads very often to self organization and to jet formation (please see for example the special 
issue of Journal of Atmospherical Science, named "Jets and Annular Structures in Geophysical Fluids" that contains the paper 
[23 1). Those jet behaviors are at the basis of midlatitude atmosphere dynamics [70] and quantifying their statistics is fundamental 
for understanding climate dynamics. A similar self-organization into jets has also been observed in two-dimensional turbulence 

0221. 

In this paper, we study the jet formation problem in the simplest possible theoretical framework: the two-dimensional 
equations for a barotropic flow with a beta effect /3. These equations, also called the barotropic quasi-geostrophic equations, 
are relevant for the understanding of large scale planetary flows [59 1. When jS = 0, they reduce to the two-dimensional Euler or 
Navier-Stokes equations. All the formal theoretical framework developed in this work could be easily extended to the equivalent 
barotropic quasi-geostrophic model (also called the Charney-Hasegawa-Mima equation), to the multi-layer quasi-geostrophic 
models or to quasi-geostrophic models for continuously stratified fluids |59|. 

The aim of this work is to consider an approach based on statistical mechanics. The equilibrium statistical mechanics of two- 
dimensional and quasi-two-dimensional turbulence is now well understood [9 62 52,63 46] : it explains self-organization and 
why zonal jets (east- west jets) are natural attractors of the dynamics. However, a drawback of the equilibrium approach is that 



Laboratoire de physique, Ecole Normale Superieure de Lyon et CNRS, 46 allee dTtalie, 69007 Lyon, France. E-mail: Freddy.Bouchet@ens-lyon.fr 
cesare.nardini@gmail.com ■ tomas.tangarife@ens-lyon.fr 



2 



Freddy Bouchet et al. 



the set of equilibrium states is huge, as it is parametrized by energy, enstrophy and all other inertial invariant of the dynamics. 
Moreover, whereas many observed jets are actually close to equilibrium states, some other jets, for instance Jupiter's ones, seem 
to be far from any equilibrium state. It is thus essential to consider a non-equilibrium statistical mechanics approach, taking 
into account forces and dissipation, in order to understand how real jets are actually selected by a non-equilibrium dynamics. 
In this work, we consider the case when the equations are forced by a stochastic force. We then use classical tools of statistical 
mechanics and field theory (stochastic averaging, projection techniques) in order to develop the kinetic theory from the original 
dynamics. 

Any known relevant kinetic approach is associated with an asymptotic expansion where a small parameter is clearly iden- 
tified. Our small parameter a [8,9] is the ratio of an inertial time scale (associated to the jet velocity and domain size) divided 
by the forcing time scale or equivalently the dissipation time scale (the spin-up or spin-down time scale, needed to reach a sta- 
tistically stationary energy balance). As discussed below, when the force is a white in time stochastic force, the average energy 
input rate is known a-priori, and the value of a can actually be expressed in terms of the control parameters . We call the limit 
a -C 1 the small force and dissipation limit. 

For small a, the phenomenology is the following. At any times the velocity field is close to a zonal jet v ~ U(y,t)e x . U (y,t)e x 
is a steady solution of the inertial dynamics. This zonal jet then evolves extremely slowly under the effect of the stochastic forces, 
of the turbulence (roughly speaking Reynolds stresses), and of the small dissipation. The non-zonal degrees of freedom have a 
turbulent motion which is strongly affected by the dominant zonal flow U(y), These turbulent fluctuations are effectively damped 
by the non-normal linear dynamics close to the zonal jets. This inviscid damping of the fluctuations is accompanied by a direct 
transfer of energy from the fluctuations to the zonal flow. The understanding and the quantification of these processes is the aim 
of this work. The final result will be a kinetic equation that describes the slow dynamics of the zonal jets, where the effects of 
the non-zonal turbulence has been integrated out. 

From a technical point of view, we start from the Fokker-Planck equation describing exactly the evolution of the Probability 
Density Function (to be more precise, Functional) (PDF) of the potential vorticity field (or vorticity field for the two-dimensional 
Navier-Stokes equations). We make an asymptotic expansion of this Fokker-Planck equation, by carefully identifying the slow 
and fast processes and the order of magnitude of all fields. At a formal level we follow an extremely classical route, described 
for instance in the theory of adiabatic averaging of stochastic processes [32], also called stochastic reduction (see also i38lfTTl 
29 1 for counterpart in mathematics). This leads to a new Fokker-Planck equation that describes the slow evolution of the zonal 
jet U ; such formal computations are tedious but involve no difficulties. 

The main theoretical challenge is to check that the asymptotic expansion is self-consistent. We have to prove that all quan- 
tities appearing in the kinetic theory remain finite, keeping in such a way the order of magnitude initially assumed. In field 
problems, like this one, this is not granted as ultraviolet divergence often occur. For example, the slow Fokker-Planck equation 
involves a non-linear force that is computed from the Ornstein-Uhlenbeck process, corresponding to the linearized dynamics 
close to U stochastically forced. This process is characterized by the two-points correlation function dynamics, which is de- 
scribed by a Lyapunov equation. We then need to prove that this Lyapunov equation has a stationary solution for a = 0, in order 
for the theory to be self-consistent. 

A large part of our work deals with the analysis of this Ornstein-Uhlenbeck process, and the Lyapunov equation, in the 
limit a — > 0. The fact that the Lyapunov equation has a stationary solution in the limit a — > is striking, as this is the inertial 
limit in which the equation for the non-zonal degrees of freedom contains a forcing term but no dissipation acts to damp them. 
The issue is quite subtle, as we prove that the vorticity-vorticity correlation function diverges point-wise, as expected based on 
enstrophy conservation. However we also prove that any integrated quantity, for instance the velocity-velocity autocorrelation 
function, reaches a stationary solution due to the effect of the zonal shear and a combination of phase mixing and global effects 
related to the non-normality of the linearized operator. As discussed thoroughly in the text, this allows to prove the convergence 
in the inertial limit of key physical quantities such as the Reynolds stress. Most of this analysis strongly relies on the asymptotic 
behavior of the linearized deterministic Euler equation, that we have studied for that purpose in a previous paper [7 ]. 

The linearized equation close to a base flow U(y) has a family of trivial zero modes that correspond to any function of y 
only. The main hypothesis of our work is that the base flows U are linearly stable (they have no exponentially growing mode), 
but also that they have no neutral modes besides the trivial zonal ones. A linear dynamics with no mode at all may seem strange 
at first sight, but this is possible for an infinite dimensional non-normal linear operator. Actually, the tendency of turbulent flows 
to produce jets that expel modes has been recognized long ago 1371 . Moreover, as discussed in [7|, the absence of non-trivial 
modes is a generic situation for zonal flows, even if it is certainly not always the case. 

We now describe the physics corresponding to the Fokker-Planck equation for the slow jets dynamics. It is equivalent to a 
stochastic dynamics for the velocity field that we call the kinetic equation. The kinetic equation has a deterministic drift F(U) 
and a stochastic force which is a Gaussian process with white in time correlation function. F(U) is the Reynolds stress that 
corresponds to the linearized dynamics close to the base flow U. If we consider the deterministic drift only, the equation is then 
related to Stochastic Structural Stability Theory (SSST) first proposed on a phenomenological basis by Farrell, Ioannou |26,25, 
[TJ, for quasi-geostrophic turbulence. It is also related to a quasi-linear approximation plus an hypothesis where zonal average 
and ensemble average are assumed to be the same, as discussed in details in section [2T4l More recently, an interpretation in 
terms of a second order closure (CE2) has also been given [49 48,69]. All these different forms of quasi-linear approximations 
have thoroughly been studied numerically, sometimes with stochastic forces and sometimes with deterministic ones [21 1. Very 
interesting empirical studies (based on numerical simulations) have been performed recently in order to study the validity of this 
type of approximation [50 48 57 69|, for the barotropic equations or for more complex dynamics. The SSST equations have also 
been used to compare theoretical prediction of the transition from a turbulence without a coherent structure to a turbulence with 
zonal jets |67 , 1 581 - Our first conclusion is that kinetic theory provides a strong support to quasi-linear types of approximations, 
in the limit of weak forces and dissipation a <C 1, in order to compute the attractors for the slow jet dynamics. 
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Beside the deterministic drift F(U), the kinetic theory also predicts the existence of a small noise. Moreover it predicts 
the Fokker-Planck equation describing the full Probability Density Functional of the zonal jet U. This was not described, even 
phenomenologically, in any previous approach. This is an essential correction in many respects. First, it allows to describe the 
Gaussian fluctuations of the zonal jet. We also note that the probability of arbitrarily large deviations from the deterministic 
attractors can be computed from this Fokker-Planck equation. For instance, we may implement large deviations theory in the 
small noise limit. This is typically the kind of result that cannot be obtained from a standard cumulant expansion or closure 
based on moments. The possibility to compute the full PDF is extremely important especially in cases where the deterministic 
part of the dynamics has more than one attractor. This case actually happens, as noticed by Ioannou and Farrell (see section[6]l. 
Then, our approach may predict the actual probability of each attractor and the transition probabilities from one attractor to the 
other. We remark anyway that, from a practical point of view, a further step forward should be done to obtain explicit results in 
this direction. 

Our work has been clearly inspired by the kinetic theory of systems with long range interactions [44 56 31, described at first 
order by Vlasov equation and at next order by Lenard-Balescu equation. We have proven in previous works that this kinetic 
theory leads to algebraic relaxations and anomalous diffusion [6, 10 73, 13]. The Euler equation is also an example of system 
with long range interaction, and there is a strong analogy between the 2D Euler and Vlasov equations. Quasilinear approximation 
for the relaxation towards equilibria of either the 2D Euler equation 1151 or the point vortex dynamics [24 16 17] have actually 
been proposed and studied in the past. 

All the above results started from deterministic dynamics, with no external forces. In order to prepare this work, and to 
extend these kinetic theories to the case with non-equilibrium stochastic forces, we have first considered the Vlasov equation 
with stochastic forces [54 55]. A kinetic equation was then obtained, similar to the one in this paper, but with much less 
technical difficulties. The reason is that the Landau damping, which is then responsible for the inviscid relaxation, can be studied 
analytically through explicit formulas at the linear level [56 44j. Non-linear Landau damping has also been recently established 
[53]. The kinetic equation for the stochastic dynamics [54 1 has the very interesting property to exhibit phase transitions and 
multistability, leading to a dynamics with random transitions from one attractor to the other [55]. We stress again that beside the 
formal structure, the main theoretical difficulty is the analysis of the Lyapunov equation which is a central object of these kinetic 
theories. In order to extend this approach to the barotropic equations, the current work discusses the first study of the Lyapunov 
equations for either the 2D Euler, the 2D Navier-Stokes or the barotropic flow equations in the inertial limit. 

The barotropic flow equations include the 2D Stochastic Navier-Stokes equations as a special case. During last decade, a 
very interesting set of mathematical works have proved results related to the existence and uniqueness of invariant measures, 
their inertial limit, their ergodicity, the validity of the law of large number 141 1140 II 3 9 II 2 81151117211 34 II 3 5 1IT21I2Q1I651I271 and of large 
deviations principles [33 68 36 j. Some of the results are summarized in a recent book [43]. In order to be applied for real physi- 
cal situations, these works should be extended in order to deal with a large scale dissipation mechanism, for instance large scale 
linear friction. We also note an interesting work considering the 2D Navier-Stokes equations forced by random vorticity patches 
1147 1 . The mathematical literature also contains a lot of interesting studies about stochastic averaging in partial differential equa- 
tions [42 1, but we do not know any example dealing with the 2D Navier-Stokes equations or the barotropic flow equations. 

In section|2j we discuss the model, the energy and enstrophy balances, non-dimensional parameters, the quasi-linear approx- 
imation and the Fokker-Planck equations for the potential vorticity PDF. We develop the formal aspects of the kinetic theory, or 
stochastic averaging, in section[3] This section ends with the derivation of the kinetic equation: the Fokker-Planck equations for 
the slow evolution of the jet ( 133b . and its corresponding stochastic dynamics, Eq. d34b or Q5b . Section|4]comes back on energy 
and enstrophy balances. The analysis of the Lyapunov equation is performed in section [5] We establish that it has a stationary 
distribution in the inertial limit, discuss the divergence of the vorticity-vorticity correlation function, the nature of its singularity 
and how they are regularized in a universal way by a small linear shear or by a small viscosity. Section [6]discusses the impor- 
tance of the stochastic part of the kinetic equation, explains how it predicts zonal jet PDF with jets arbitrarily far from Gaussian 
fluctuations, and stress the existence of cases with multiple attractors, phase transitions and bistability. We discuss open issues 
and perspectives in section|7] The paper contains also three appendices reporting the most technical details of sections[3]and|5] 



2 Quasi-Geostrophic dynamics and Fokker-Planck equation 

2. 1 2D barotropic flow upon a topography 

We are interested in the non-equilibrium dynamics associated to the 2D motion of a barotropic flow, with topography h, on a 
periodic domain @ = [0,2kI x L) x [0,2kL) with aspect ratio l/l x - 



where CO, q, v and iff are respectively the vorticity, potential vorticity, the non-divergent velocity, and the stream function. 
In the equations above, r = (x,y) are space vectors, x is the zonal coordinate and y the meridional one, A is a Rayleigh (or 
Ekman) friction coefficient and v n j is the hyper- viscosity coefficient (or viscosity for n = 1). All these fields are periodic 
f(x + 2Kl x L,y) = f(x,y) and / (x,y + 2kL) = f(x,y). We have introduced a forcing term rj, assumed to be a white in time 
Gaussian noise with autocorrelation function E [Tj(ri,fi)T}(r2,?2)] = C(ri — T2)5(fi — tz), where C is an even positive definite 
function, periodic with respect to x and y. As discussed below, a is the average energy input rate. When h = 0, the 2D barotropic 




— +v Vq = -Xto-v, ud (-A) n co + \fdr\, 
at 

y = e,xV^, co = A\j/, q = co + h(y), 



(1) 
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flow equations are the 2D-Navier-Stokes equations. Here we assume that the noise autocorrelation function C is translationally 
invariant in both direction. The fact that C is zonally invariant is important for some of the computations in this paper. However 
the hypothesis that C is meridionally invariant is not important and generalization to non-meridionally invariant forcing would 
be straightforward. 

Dynamical invariants of perfect barotropic flows Equations (TJl with X = a = v n j = describe a perfect barotropic flow. The 
equations are then Hamiltonian and they conserve the energy 

S[q] = \f drv 2 = -\ [ dr (q-h) ¥ , (2) 
2 J 9 2 J s< 

and the Casimir functionals 

%\q] = / dr s(q), (3) 
for any sufficiently regular function s. The potential enstrophy 

¥ 2 [q\ = ± { Arq 2 
2 .19 

is one of the invariants. When h = 0, the perfect barotropic flow equations clearly reduce to the 2D Euler equations. 



Averaged energy input rate and non-dimensional equations Because the force is a white in time Gaussian process, we can com- 
pute a-priori the average, with respect to noise realizations, of the input rate for quadratic invariants. Without loss of generality, 
we assume that 

-2n 2 l x L 2 {A- l C) (0) = 1, 

where A ~ 1 denotes the inverse Laplacian operator; indeed, multiplying C by an arbitrary positive constant amounts at renormal- 
izing a. Then, with the above choice, the average energy input rate is a and the average energy input rate by unit of mass is 
e = a j\n 2 l x L 2 . Moreover, the average potential enstrophy input rate is given by 

2ti 2 I x L 2 C{0)o. 



We consider the energy balance for equation ([T}, with E = E [t£%]]: 




dE 

— = -2XE-v n , d H n + o, (4) 
dt 

where H„ = — E \j/(—A)" co\. For most of the turbulent flows we are interested in, the ratio 2XE/v„ jH n will be extremely 
large (viscosity is negligible for energy dissipation). Then, in a statistically stationary regime, the approximate average energy is 
E ~ o/2X. We perform a transformation to non-dimensional variables such that in the new units the domain is $> = [0,2tzI x ) x 
[0, 2k) and the approximate average energy is 1. This is done introducing a non-dimensional time t' = tjx and a non-dimensional 
spatial variable r' = r/L with T = L 2 \/2X /a. The non-dimensional physical variables are q 1 = Tq, v' = Tv/L, h' = %h, and the 
non-dimensional parameters are defined by 



a = At = L 



v„ = v n ^xjL 2n = v n/ 2 \j2Xja jL 2n 2 . We consider a rescaled stochastic Gaussian field r)' with E {v\,t'^)i]{v' 1 ,t' 7 )\ = C'(r[ — 
r' 2 )5(r' 1 - 4) with C'(r') = L 4 C(r). Performing the non-dimensionalization procedure explained above, the barotropic equations 
read 

{^+v Vq= -aco- v„(-A)" w + V2arj : 
at (5) 
v = e ; xVij/, a = Ay, q=a + h(y), 

where, for easiness in the notations, we drop here and in the following the primes. We note that in non-dimensional units, a 
represents an inverse Reynolds number based on the large scale dissipation of energy and v„ is an inverse Reynolds number 
based on the viscosity or hyper-viscosity term that acts predominantly at small scales. The non dimensional energy balance is 

dE 

— = -2aE - v n H„ + 2a. (6) 
dt 
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2.2 2D barotropic flow with a beta-effect 

For a 2D flow on a rotating sphere, the Coriolis parameter depends on the latitude. For mid-latitudes, such a dependence can be 
approached by a linear function, the so-called beta effect: this corresponds to consider a linear topography function h(y) = jS^y. 
In a periodic geometry, such a function is meaningless, and so is the associated potential vorticity q = co + h. However, consistent 
equations can be written for the vorticity CO and the velocity: 

{^+v- Vw+pvy = -aco-v„(-A)" co + Vlari, 
dt (7) 
v = e z xVi^, co = A\j/ 

Observe that, in this case, we have directly written the equations in a non-dimensional form; the parameter jS as a function of 
dimensional quantities is given by 

2T„ / L x 1 



P=L\ — p d 



where /3j is the dimensional beta effect, and we have defined a Rhines scale Ln = a/L/t/^ based on the large scale velocity 
U = L/t. 

The perfect barotropic equations (a = v„ = 0), in doubly periodic geometry, with beta effect (jS ^ 0), conserve only two 
quantities: the kinetic energy 



<f[fl>] = -~ / draw, (8) 

2 JS 



and the enstrophy 



3f\co] = - [ dr co 2 . (9) 

2 J 9 

For jS = 0, the beta-plane equations reduce to the 2D-Navier Stokes equations discussed in previous section. 

For sake of simplicity, in the following, we consider the case of a viscous dissipation n = 1, and we denote the viscosity 
v = Vi. One of the consequences of our theoretical work is that, in the limit a< 1 and v„ <C cc, the results will be independent 
on v„ for generic cases, and at leading order in the small parameters. 



2.3 Decomposition into zonal and non-zonal flow 

The physical phenomenon we are interested in is the formation of large scales structures (jets and vortices). Such large scale 
features are slowly dissipated, mainly due to the friction a. This dissipation is balanced by Reynolds stresses due to the transfer 
of energy from the forcing scale until the scale of these structures. This phenomenology is commonly observed in planetary 
atmospheres (Earth, Jupiter, Saturn) and in numerical simulations. For the barotropic equations (|5} and (|7}, the regime corre- 
sponding to this phenomenology is observed when acl (time scale for forcing and dissipation 1 /A much larger than time 
scale for the inertial dynamics t) and when 2XE /v n jH n 2> 1 (turbulent regime). For this reason, we study in the following the 
limit v„ — > 0, a — > (a <C 1 and v„ <C a). 

For sake of simplicity, we consider below the case when the zonal symmetry (invariance by translation along x) is not broken. 
Then the large scale structure will be a zonal jet characterized by either a zonal velocity field v z = U(y)e x or its corresponding 
zonal potential vorticity q z (y) = —U'(y) + h(y). For reasons that will become clear in the following discussion (we will explain 
that this is a natural hypothesis and prove that it is self-consistent in the limit a <C 1), the perturbation to this zonal velocity field 
is of order y/a. 

Defining the zonal projection (.) as 

the zonal part of the potential vorticity will be denoted by q, = (q) ; the rescaled non-zonal part of the flow q m = co m is then 
defined through the decomposition 

q{r)=q z (y) + \/aco m {r). 

The zonal and non-zonal velocities are then defined through U'(y) = —q z {y) +h(y), the periodicity condition, and 

v(r) = U(y)e x + ^a\ m (r). 
We now project the barotropic equation (|5} into zonal 

dq? d I f v j \ d 2 C0 7 r— 

and non-zonal part 



^-+Lu [co m ] = \flr\ m - x/a\, n .Vco m + Va {\ m .Vco,„) , (1 1) 
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where r\, = (r\) is a Gaussian field with correlation function E [T/ z (yi,fi)7J z (y2,f2)] = ^z(yi ~y2)8(h —t%) with C z = (C), 
Tim = tj — (tj) is a Gaussian field with correlation function E [T/ m (ri,fi)T/ m (r2,f2)] — C m (r\ —i2)8(t\ — tz) with C m = C — (C), 
and where 

L u [(o m ]=U(y)-^-+q',{y)-^- + aco m -vAco m , <D m =Ay m . (12) 
ox ox 

We observe that, in the previous equations, —jf 1 +Lu [co m ] = is the deterministic dynamics linearized close to the zonal base 
flow U, whose corresponding potential vorticity is q z . In the following we will also consider the operator 

rO r,-. 1 Til \ a ®m . I, x^Vm ,,,,, 

L u [co m \ = U{y)^—+q z (y)^—, (13) 
ox ox 

which is the operator for the linearized inertial dynamics (with no dissipation) close to the base flow U. In all the following, 
we assume that the base flow U(y) is linearly stable (the operator has no unstable mode). We remark that the action of Lfj 
on zonal functions f{y) is trivial: any zonal function is a neutral mode of L° . We will thus consider the operator acting 
on non-zonal functions only (functions which have a zero zonal average). We assume that L° restricted to non-zonal functions 
has no normal mode at all (which is possible as L^, is a non-normal operator). This hypothesis will be important for the results 
presented in section [5] This hypothesis may seem restrictive, but as explained in 1 7 1 this is a generic case, probably the most 
common one. 



The equation for the zonal potential vorticity evolution ( 110b can readily be integrated in order to get an equation for the zonal 
flow evolution 

du I (>') \ r , d2 U 

-fa= a ym®mj -ccU + + V2a-qu, 

where r\u is a Gaussian field with correlation function ^[T]u{y\ih)T]u{y2ih)\ = Cu(y\ — yz)8(t\ — tz) with ^r(y\ — yi) = 

Qfji -yz)- 

We see that the zonal potential vorticity is coupled to the non-zonal one through the zonal average of the advection term. 
If our rescaling of the equations is correct, we clearly see that the natural time scale for the evolution of the zonal flow is I /a. 
By contrast the natural time scale for the evolution of the non-zonal perturbation is one. These remarks show that in the limit 
a <C 1 , we have a time scale separation between the slow zonal evolution and a rapid non-zonal evolution. Our aim is to use 
this remark in order to i) describe precisely the stochastic behavior of the Reynolds stress in this limit by integrating out the 
non-zonal turbulence, ii) prove that our rescaling of the equations and the time scale separation hypothesis are self-consistent. 

The term y / av m .V co m — y/a{\ m .V(O m ) describes the interactions between non-zonal degrees of freedom (sometimes called 
eddy-eddy interactions). If our rescaling is correct, these terms should be negligible at leading order. Neglecting these terms 
leads to the so called quasi-linear dynamics, which is described in next section. 



2.4 The quasi-linear dynamics 

By neglecting the interactions between non-zonal degrees of freedom in (lilt , one obtains the quasi-linear approximation of the 
barotropic equations: 

dq z d I (yi \ d 2 (£>, r— 

St dy\ I dy z (14) 

d(O m r- 

-^-+Lu [com] =V2r\ m 

In two-dimensional and geostrophic turbulence, the fluctuations around a given mean flow are often weak, so this approximation 
is natural in this context. If our rescaling is relevant, this corresponds to the limit acl. 

The approximation leading to the quasi-linear dynamics amounts at suppressing some of the triads interactions. As a con- 
sequence, the inertial quasi-linear dynamics has the same quadratic invariants as the initial barotropic equations: the energy and 
potential enstrophy. 

As discussed in the previous paragraphs, and as can be seen in J 14b . in the regime a <C 1, it seems natural to assume that there 
is a separation between the time scale for the fluctuation dynamics (second equation of l!14t). which is of order 1, and the time 
scale for the zonal flow (first equation of (I14t ). which is of order 1 /a. At leading order, the evolution of co m can be considered 
with the zonal velocity profile held fixed. The second equation of d 1 4b . with a fixed velocity profile, is then linear. Thus, denoting 
by E the average over the realization of the noise r\ m for fixed U, the distribution of co m is completely characterized by the two- 
points correlation function g(ri,r2,f) = E [co m (r\ ,f)<» m (r2,f)]. The evolution of g is given by the so-called Lyapunov equation, 
which is obtained by applying the Ito formula (with U fixed): 

where L,ff (resp. iff) is the linearized operator Lu J 12b acting on the variable ri (resp. Vi). 

If we assume that this Gaussian process has a stationary distribution, then considering again the time scale separation a< 1, 
one is led to the natural hypothesis that at leading order 

dq z d I r v i \ d 2 CO, i — 

= -a-^-Eu {vm'cOmj - aco z + v-^f- + V2arj z , (15) 
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where Ej/ (vm" 1 03m) is the average over the stationary Gaussian process corresponding to the dynamics given by the second 

equation of J 14b at fixed U. In our problem, invariant by translation along the zonal direction, we note that — Ej/ (v^ton,^ is 
the divergence of a Reynolds stress [60 1. We call this equation the kinetic equation, because it is the analogous of the kinetic 
equation obtained in a deterministic context in the kinetic theory of plasma physics [56 ..441, systems with long-range interactions 
[5,6 10] and the 2D Euler equations 1150161 . The time scale separation hypothesis, leading to consider the asymptotic solution 
of the Lyapunov equation is referred to as the Bogolyubov hypothesis in classical kinetic theory 15611441 . 
Another model related to the quasi-linear dynamics is 



dq 7 d [/ M \1 d^W, i — 

= -a-^E,„ l^vJ^co^J -aco z + v-^f + V2arj z 

d 4 + (Ly+L$) g = 2C m 



u 7 m — 
Y + v 2ar\ z 

(16) 

dt 

where E m is the average corresponding to the Gaussian process with two-point correlation function g(r, r', t) = E„, [co„, (r , t ) co m (r', ?)]. 
By contrast to the kinetic equation ( 115b . in equation l l!6[ i g and U (or q z ) evolve simultaneously. The Reynolds stress divergence 
— E„, ^Vm^O) m ^j , or the gradient of Reynolds stress divergence — J^E m j^v„'a) m ^ j , can be computed as a linear transform of 
g. This model, first proposed by Farrell and Ioannou [26 25 fl, is referred to as Structural Stochastic Stability Theory (SSST). 
Farrell and Ioannou have also used this model in a deterministic context [21 1, where the correlation C m is then added phenomeno- 
logically in order to model the effect of small scale turbulence (the effect of the neglected non linear terms in the quasi-linear 
dynamics). This model, or an analogous one in a deterministic context, is also referred to as CE2 (Cumulant expansion of order 
2) B0ll49ll4in . 

As previously observed, the second equation in ( 1161 can be deduced from the second equation of J141 as an average over 
the realizations of the non-zonal part of the noise T} m , with held fixed and independent on U (neither of which are true). It can 
be considered as an approximation of the quasilinear dynamics, where the instantaneous value of the non-linear term / v~m (0„ 



is replaced by its ensemble average, with U held fixed. We see no way how this could be relevant without a good time scale 
separation a <C 1. With a good time scale separation, both the quasilinear dynamics ( 1141 and SSST-CE2 d 16b are likely to be 
described at leading order by the kinetic equation (1151 . We see no reason how SSST dynamics could be relevant to describe the 
quasilinear dynamics beyond the limit where it is an approximation of the kinetic equation. However, from a practical point of 
view, SSST-CE2 is extremely interesting as it provides a closed dynamical system which may be extremely easily computed 
numerically (by comparison to both direct numerical simulations, the quasilinear dynamics or the kinetic equation) in order 
to obtain the average zonal velocity and its Gaussian corrections. It is probably the best tool to understand qualitatively the 
dynamics of barotropic flows in the limit of small forces and dissipations. 



2.5 Fokker-Planck equation 

We will use the remark that we have a time scale separation between zonal and non-zonal degrees of freedom in order to average 
out the effect of the non-zonal turbulence. This corresponds to treat the zonal degrees of freedom adiabatically. This kind of 
problems are described in the theoretical physics literature as adiabatic elimination of fast variables 1321 or stochastic averaging 
in the mathematical literature. Our aim is to perform the stochastic averaging of the barotropic flow equation and to find the 
equation that describes the slow evolution of zonal flows. 

There are several ways to perform this task, for instance using path integral formalism, working directly with the stochastic 
equations, using Mori-Swanzig or projection formalisms, and so on. A typical way in turbulence is to write the hierarchy of 
equations for the moments of the hydrodynamic variables. As discussed in the introduction and in section [6] this approach may 
lead to interesting results sometimes, but it can sometimes be misleading, and when it works it can describe quasi-gaussian 
fluctuations at best. The formalism we find the simplest, the more convenient, and the more amenable to a clear mathematical 
study is the one based on the (functional) Fokker-Planck equations. 

At a formal level, we will perform adiabatic reduction of the Fokker-Planck equation using the classical approach, as de- 
scribed for instance in Gardiner textbook [32 1. Whereas Gardiner treats examples with a finite number of degrees of freedom, 
we are concerned in this paper with a field problem. At a formal level, this does not make much difference and the formalism 
can be directly generalized to this case. At a deeper level this however hides some mathematical difficulties, some of them being 
related to ultraviolet divergences. We will show that such ultraviolet divergence do not occur, at least for the quantities arising 
at leading order. 

The starting point of our analysis is to write the Fokker-Planck equation associated with Eq. d 1 Ob and i ll lb . We consider the 
time dependent probability distribution function P, [q z , co m \ for the potential vorticity field and, for easiness in the notations, we 
omit time t in the following. The distribution P [q Z7 co m \ is a functional of the two fields q z and co,„ and is a formal generalization 
of the probability distribution function for variables in finite dimensional spaces. From the stochastic dynamics d 101 lib , using 
standard techniques (functional Ito calculus), one can prove that P solves the Fokker-Planck equation 

dP — 

-zp = ^f Q P + s/aJf„P + a^-P, (17) 
at 

where -g-^yj are functional derivatives with respect to the field q at point y, where the leading order term is 



Sf P= / dn-= — 

<5co m (ri 



dP 

L u [a m ](r x )P+ / dr 2 C m (ri-r 2 )^ -- 

0(D m {r 2 ) 



(18) 
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the zonal part of the perturbation is 

d (vlZ ] (O m ^{yi) + CO z {yi)-^Aco z {yi))p + J dy 2 C z (yi-y 2 ) 8 ' 



&,P = I dvi ' 



dy \ I a J J oq z (y 2 



(19) 



and the nonlinear part of the perturbation is 



dri ^ — -— [(v m .V© m (n) - (v„,.Vffl m (n)))P] . (20) 
oco m (ri) 

We observe that J£ n is a linear operator which describes the effect of the non-linear terms in the fluid mechanics equations. 

The operator J£q is the Fokker-Planck operator that corresponds to the linearized dynamics close to the zonal flow U, forced 
by a Gaussian noise, white in time and with spatial correlations C m . This Fokker-Planck operator acts on the non-zonal variables 
only and depends parametrically on U. This is in accordance with the fact that on time scales of order 1, the zonal flow does 
not evolve and only the non-zonal degrees of freedom evolve significantly. It should also be remarked that this term contains 
dissipation terms of order a and V. These dissipation terms could have equivalently been included in Jf z , but we keep them on 
Jz?o for latter convenience. However it will be an important part of our work to prove that in the limit v< a< 1, at leading 
order, the operator Jzfo with or without these dissipation terms will have the same stationary distributions. This is a crucial point 
that will be made clear in section [5] 

At order y/a, the nonlinear part of the perturbation Jzf„ contains the non-linear interactions between non-zonal degrees 
of freedom. We will prove that at leading order Jzf„ has no effect, justifying the quasilinear approach. This is non-trivial as Jzf„ 
formally arise at an order in \/a lower than J5? z . At order a, the zonal part of the perturbation Jzf z contains the terms that describe 
the coupling between the zonal and non-zonal flow, the dynamics due to friction acting on zonal scales and the zonal part of the 
stochastic forces. 



3 Stochastic averaging and the Fokker-Planck equation for the slow evolution of zonal velocity profiles 

In this section, we formally perform the perturbative expansion in power of yfa, for the Fokker-Planck equation dl7t . The aim of 
the computation is to obtain a new equation describing the slow zonal part only. We will compute explicitly only the leading order 
relevant terms. This formal computation will make sense only if the Gaussian process corresponding to a stochastically forced 
linearized equation has a stationary solution. That last point is thus the real physical issue and the most important mathematical 
one. We consider this question in section|5] 



3.1 Stationary distribution of the fast non-zonal variables 
At leading order when a is small, we have 

%=^R (21) 

Let us first consider the special case when the zonal flow is a deterministic one: P[q,co m ] = 8(q — q z )Q(co m ). Then, ( 121b is 
a Fokker-Planck equation corresponding to the dynamics of the non-zonal degrees of freedom only. It is the Fokker-Planck 
equation associated to the barotropic equations linearized around a fixed base flow with zonal velocity U and forced by a 
Gaussian noise delta correlated in time and with spatial correlation function C,„: 

^-+L u [a m ] = V2r\ m , E[Tj m (ri,ri)7] w (r 2 ,fc)] = C m (n -r 2 )S(fi -t 2 ) (22) 
at 

When U is held fixed, this is a linear stochastic Gaussian process, or Ornstein-Uhlenbeck process. Thus, it is completely char- 
acterized by the two-points correlation function g(ri,r 2 ,f) = IE [co m (ri ,f)a>„, (r2,f)] (the average E[a>,„] is equal to zero, here the 
average E refers to an average over the realization of the noise r\ m , for fixed U). The evolution of g is given by the so-called 
Lyapunov equation, which is obtained by applying the Ito formula to l!22t : 

||+L* / 1) g + 4 2) g = 2C m , (23) 

where Ly (resp. Ly ) is the linearized operator Ly d 12b acting on the variable ri (resp. r 2 ). 

For fixed a and v, the linear operator Ly fl!2b is dissipative. If the linearized dynamics is stable, then the Ornstein-Uhlenbeck 
will have a stationary distribution. However, as we are computing an asymptotic expansion for small a and v, we need to consider 
the limit V — > and a — > of this stationary distribution. In order to do that we will consider the Lyapunov equation 



df - +L y W g + Lf ) g = 2C m , (24) 



where Ly is the inertial linear dynamics d!3t . In this section, we assume that the corresponding Gaussian process has a stationary 
distribution. This assumption may seem paradoxical as there is no dissipation in L°; however, it will be proved in section|5]that 
this is correct. We denote by g°° [q z ] the stationary two-points correlation function, g°°[<fe](ri,r 2 ) = lim,^ 00 g(ri ,r 2 ,f), and by 
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(s°°) ' [iz] (1*1 , 1*2) its invers^]. Then, the stationary distribution of the linear equation (122b close to the base flow with potential 
vorticity q z is 



G[q z ,eo m ]= -j== exp\--J dridr 2 (0 m (r 1 )(g 00 ) [q z ] (ri,r 2 )co m (r 2 )J , (25) 

where the normalization constant depends on detg., the determinant of 2ng°° [q z ]. The stationary solution to d21b with any initial 
condition P[q, C0 m ] = 8(q — q z ) Q(co m ) is thus given by 8 (q— q z )G[q z , C0 m ]. 

For any q z , the average of any observable A [»„,] over the stationary Gaussian distribution will be denoted 



E v [A] = J & [co m ] A [co m ] G [q z , w m ] 
and, for any two observables A [a>,„] and B [»„,], the correlation of A at time t with B et time zero will be denoted 

E v [A(r)S(0)] = 1 9 [co m ] A [©,„] e'^» [B [co m ] G [q z , co m }} . (26) 

The covariance of A at time t with B at time zero is denoted by 

Eu [[A(f)£(0)]] = E v {(A(t) - Ev [A]) (5(0) - E v [5])] . (27) 

As already pointed out, we assume that for all q, of interest the Gaussian process corresponding to the inertial linear dynamics 
d 1 3 b has a stationary distribution. Then, the operator exp (f-Sfrj) has a limit for time going to infinity. We can then define 9* the 
projector on the asymptotic solutions of this equation 

,0 s [P] = Kmexp(f-2o)^ 

for any P. We observe that, because J£q does not act on zonal degrees of freedom, [8 (q z — q) Q((O m )\ = 8 (q z — q)G [qo, co m \ 
for any Q. Moreover, using the linearity of we conclude that for any P 

9>\P\ = G \q z , co m \ J $ [cow,] P [q z , a>m] ■ (28) 

The above equation expresses that for any q z , fast evolving non-zonal degrees of freedom relax to the stationary distribution 
G. With this last relation, it is easily checked that 9 is a projector: 3? 1 = Moreover, because commutes with ^fo and 
projects on its stationary distribution, we have 

^^0 = ^0^ = 0. 



3.2 Adiabatic elimination of the fast variables 

We decompose any PDF P into its slow and fast components: 

P = P s +P f 

with P s = 9>P and P f = (1 - &>)P. 
Using ( 128b . we obtain 

Ps[qz,tOm] = R[qz\G[q z ,(O m ], 

with R [q z ] = J & [co m ] P s \q z , CO m ]. R is the marginal distribution of P s on the zonal modes and describes the statistics of the zonal 
flow; P s describes the statistics of the zonal flow assuming that for any q z , the fast non-zonal degrees of freedom instantaneously 
relax to their stationary Gaussian distribution (I25t . We thus expect R and P s to evolve slowly, and Pf contains the small correc- 
tions to P s . The goal of the stochastic reduction presented here is to get a closed equation for the evolution of R, valid for small 
a. 

Applying the projector operator 9* to Eq. d 1 7b and using 3?J£q = 0, we get the evolution equation for P s 

DP — 

-j- = 9 + VoL&h) (Ps + Pf) ■ 

As expected, this equation evolves on a time scale much larger than 1, and this is due to the relation ^«5?o = 0, which is 
essentially the definition of 9*. 

Using ( 128b and ( 120b we remark that 5?I£ n involves the integral over 3$[co m ] of a divergence with respect to co m . As a 
consequence, 

3?% n = 0. (29) 

An important consequence of this relation is that the first-order non-linear correction to the quasi-linear dynamics is exactly 
zero: 

dP 

= a @>££ z (P s + P f ) . (30) 

1 Here, inverse is understood in the linear operator sense: J dr'g°°[go]( r l ,r') (g°°) 1 [^o](r',i"2) = S(l\ — 12). 
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Applying now the operator (1 - S 3 ) to Eq. (17]|, using (1 - ^)Jzfo = Sfo(l — and (|29]l, we obtain 

t)Pf 

-jj- = J? Pf + (VZJ?„ + ct( 1 - £»)J2J) (P s + P f ) . 



(31) 



Our goal is to solve ( 131b and to inject the solution into ( DOb . which will become a closed equation for the slowly evolving PDF 
P s . We solve equations ( 130131b using Laplace transform. The Laplace transform of a function of time / is defined by 



/(*) = / dre-*/(0, 
Jo 



where the real part of s is sufficiently large for the integral to converge. Using (d,f)(s) = sf(s) — /(0) and the fact that the 
operators don't depend explicitly on time, equations ( 130b and ( 131b become 



sP s (s) = a&££ z (fi+Pf) +P S (0) 

sP f {s) = Jf Pf + (o(l - &)Sf t + VaSf„) (P s + Pf)+P f {0). 



For simplicity, we assume that the initial condition satisfies P/{0) = 0. In most classical cases, this assumption is not a 
restriction, as relaxation towards such a distribution is exponential and fast. For our case, this may be trickier because as discussed 
in section [5] the relaxation will be algebraic for time scales much smaller than I /a. However, we do not discuss this point in 
detail before the conclusion (section [7}. 

Denoting by -g> the inverse of a generic operator jSf , the second equation can be solved as 



Pf ^ v "*-^f -\/a-%-ce(l 
and this solution can be injected into the first equation: 

5e n + s[a(\-&>)^ 



sP s (s) 



a^^ + a^ 2 ^^- 



Va^-afl-Vr 



Ps + Ps{0). 



At this stage we made no approximation and this last result is exact. 

We look at an expansion in powers of ^/a for the above equation. At order 4 in ^/a, we get 



sP s (s) = \ a^Jzf, + o? l2 3?5e z — l -T2r5? n 



P.+P. 



(0) + ^(a 5 / 2 ) 



We use that s Jg> is the Laplace transform of exp(?Jzfo), that I S J^> ) is the Laplace transform of fexp(r«£o), and that the 
inverse Laplace transform of a product is a convolution to conclude that 

dPs 



=a&S£ t P s + a 3/2 ^^ j dt' y^% n P s (t - f')] 



+ a 2 St>^ z £ dt' [( 1 - &>)% z + 1 '^ 2 ] P s (f - 1' ) + (a 5/2 ^j 



We observe now that the evolution equation for P s contains memory terms. However, in the limit a <C 1, P s evolves very slowly. 
Thus, if we make a Markovianization of the time integrals replacing P s (t — t) by P s (t), we make an error of order ~& (t) which 
is of order a. For example, the first integral in the last equation is 



' dt'e'''^5?„P s (t) + 0(a). 



The evolution equation for P s is then 
dP s 



dt 



= < a 



^^ z + a? l2 ^^ z £ ' dt' + a 2 &£? z £ dt' e''^ [(1 - &>)Sf z + t'^ 2 } j P s + G (a 5/2 ) . (32) 



We now have a closed differential equation for the slow probability distribution function P s . 
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3.3 The Fokker-Planck equation for the slow evolution of the zonal flow 

The explicit computation of each term involved in ( 132b . which is reported in the appendix [A] leads to the final Fokker-Planck 
equation for the zonal jets: 



3s = /a„ 5 



dz J 8q z (yi) 



dFi [U] , , V d 2 w- 



Rq : ! j dy 1 j^— ) {C R {y u y 2 )[q z ]R[q z ])\ - (33) 



dy\ ' a dy 

which evolves over the time scale T = at, with the drift term 

F X [U]=F[U] + aJ°°dt't'E v [(v^flfe) (yi,t')M[q z ,w m ]{Q)\ , 

with M given by Eq. l |92| | in appendixfAj with 

F[U]=E u \ l ^ ) co m )'j (yi), 

and with the diffusion coefficient 

C R (y u y 1 )[q z }=C z (y 1 -y 2 ) + a£dt'j^ r E u [[(v^ahn) (yi,0 (v£ y W) feO)]] . 

From its definition, we see that F [U] is the opposite of the Reynolds stress divergence computed from the quasi-linear approxi- 
mation. 

The Fokker-Planck equation ( 133b is equivalent to a non-linear stochastic partial differential equation for the potential vorticity 

where £ is a white in time Gaussian noise with spatial correlation Cr. As Cr depends itself on the velocity field U, this is a 
non-linear noise. 

At first order in a, we recover the deterministic kinetic equation (115) discussed in section l2~4l The main result here is that the 
first contribution of the non-linear operator, the order 0(a 3 / 2 ) in 02t . is exactly zero. We could have applied the same stochastic 
reduction techniques to the quasi-linear dynamics (114) and we would have then obtained the same deterministic kinetic equation 
at leading order, as expected. We thus have shown that at order 3 in \/a, the quasi-linear approximation and the full non-linear 
dynamics give the same results for the zonal flow statistics. 

At next order, we see a correction to the drift F[U] due to the non-linear interactions. At this order, the quasilinear dynamics 
and non-linear dynamics differ. We also see the appearance of the noise term, which has a qualitatively different effect than 
the drift term. We note that at this order, we would have obtained the same non-linear noise Cr from the quasilinear dynamics. 
Pushing our computations to next orders, we would obtain higher order corrections to the drift and noise terms, for instance due 
to eddy-eddy (non-zonal) non-linear interactions, but also corrections through fourth order operators. 

We note that we can also write directly an equation equivalent to 04b for the velocity profile 

^=F l [U] + i;, (35) 



where is a white in time Gaussian noise with spacial correlation 

% [§(yi,ri)SC».fe)]= (c u (y l -y 2 ) + aJ~dt'E u [[(v£W) M (v£W) C»,0)]]) 8{ h - h ) 



(36) 

Whereas at a formal level, correction to Fq and the noise term appear at the same order, their qualitative effect is quite different. 
For instance if one is interested in large deviations from the most probable states, correction of order a to Fq will still be 
vanishingly small, whereas the effect of the noise will be essential. At leading order, the large fluctuations will be given by 

§|=F[t/] + |, (37) 

Equation ( 137b then appears to be the minimal model in order to describe the evolution of zonal jet in the limit of weak forcing 
and dissipation. We will comment further on this issue in section|6l dealing with bistability of zonal jets and phase transitions. 



4 Energy and enstrophy balances 

We discuss here the energy balance in the inertial limit a <C 1 and the consistency of the stochastic reduction at the level of the 
energy. It is thus essential to distinguish the different ways to define the averages, for the original stochastic equation ^ or for 
the zonal Fokker-Planck equation ( 133b . We recall that E[-] = / £>[q z ]3>[(0 m ] ■ P[q z , C0,„] denotes the average with respect to the full 
PDF P, whose evolution is given by the Fokker-Planck equation l !17b . or equivalently, the average over realizations of the noise 
rj in equation J 1011 1 b ; Ey [■] = / &[(O m ] ■ G[q z , co m ] denotes the average with respect to the stationary Gaussian distribution of the 
non-zonal fluctuations G, defined in ( 125b and Er[-] = J &[q z \ -R[q z \ denotes the average with respect to the slowly evolving zonal 
PDF R whose evolution is given by the zonal Fokker-Planck equation 03b . Er is equivalently an average over the realizations 
of the noise £ in equation ( 134b . 

As discussed in the presentation of the barotropic equations, we are interested in the regime where the dissipation of energy 
is dominated by the one due to the large-scale linear friction: 2XE /v n jH n 2> 1 in equation As a consequence, we will 
consider in this section only the case of zero viscosity, v = 0. 
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4. 1 Energy balance for the barotropic equations 

The total energy balance of the non-dimensional barotropic equations (|5} given by Eq. (f6]l is reported here for convenience: 

dE 

— = -2aE + 2a (38) 
at 

where we recall that E = E [<o°] . With the orthogonal decomposition into zonal and non-zonal degrees of freedom, we have a 
natural decomposition of the energy contained in zonal and non-zonal degrees of freedom: E = E Z + aE m , with E z = E [5 f@ U 2 ] 
andE m = E[^j^\l]. 

Zonal energy balance From the definition of E z , either by direct computation from the Fokker-Planck equation l !17b or from Eq. 
d 1 Ob applying the Ito formula, we have 

^ = a2nl x JdyE [(v$ co„^ (y)U(y)j - 2aE z + 2aa z (39) 

with a z = -2k 2 l x (A -1 Q) (0) the rate of energy injected by the forcing directly into the zonal degrees of freedom. In addition 
with the expected energy dissipation and direct energy injection by a z , the first term on the right hand side describes the energy 
production due to the non-zonal fluctuations. 



Non-zonal energy balance Using E = E z + ccE m and equations d38b . 039b . we obtain 

^ = -2nl x J dyE [(v^ v) a) m ) (y)U(y)\ - 2aE m + 2a m , (40) 

where a m = 1 — O z is the rate of energy injection by the forcing on non-zonal degrees of freedom. Clearly, Eq. ( 139) and ( 140b are 
exact and fully equivalent to the energy balance l l3"8l >. 



4.2 Energy and enstrophy balance for the kinetic equation 



We now show that the energy balance of the kinetic equation ( 134b is consistent with the exact energy balances written above, 
at leading order in a. We denote by E z = Er [j f@ U 2 ] the average zonal energy for the kinetic equation and by a£ m = 
aEu [jj^v^] the energy contained in non-zonal degrees of freedom. We start from the energy balance for zonal degrees 
of freedom. Again, working at the level of the zonal Fokker-Planck d33b equation or of the kinetic equation d34b give the same 
result: 

rip r r 

' " ■ ~ J dyF l [U}(y)U(y) 



dx 



-2nl x E R 



-2E z -2tzI x E r 



dy (A-'Cr) (0) 



To approximate the above equation at leading order in a, we observe that F% [U] = F[U] + 0(a) = Ej/ 
and Cr = C z + 0(a). Then, the stationary energy balance for zonal degrees of freedom reduces to 



(yi) + 0(a), 



dE, 
dx 



-2nl x E R 



dvE, 



V 



■2E z + 2a z + 0{a) 



(41) 



In the limit a — > 0, the full PDF P[q z ,co„,\ is given by the slowly evolving part P s [q z ,co m ] 
energy transferred from the fluctuations to the zonal flow becomes 



G[q z , (O m ]R[q z ]. Then, the rate of 



E« 



J AyE v [(v^'ffl,,,)] U +0{a) = J dyE [(v^a^U 



and equation d4 1 b reduces to the energy balance for zonal degrees of freedom ( 139b at leading order in a. This proves the 
consistency of the zonal energy balance for the kinetic equation and the barotropic equations: E z = E- + 0(a). 



4.2.1 Non-zonal energy balance and total transfer of energy to the zonal flow 

Let us now consider the energy balance for aE m . Applying Ito's formula to J22b , we get, for any U, 

= -2tiI x J dyEu \ i ^; ) co m )(y)~\u(y)-2aE m + 2o m , (42) 

with a m = —2n 2 l x (4 _1 C,„) (0). Equation (142b does not contain time evolution which is consistent with the definition of E m in 
the kinetic equation through a stationary average. In the limit a <C 1, in agreement with the scaling of the variables we expect 
the energy of the non-zonal degrees of freedom to be of order a, or E m = 0(1). This is essential for the consistency of the 
asymptotic expansion and will be proved in section |5] 

From E m = 0(1) and equation d42b . the energy dissipated in the non-zonal fluctuations per unit time is negligible, so the 
stationary energy balance for the non-zonal fluctuations gives 

nl x JdyE [(v^co,,,) u] = M x J dyE v [(v&W) (y)} U(y) + 0(a) = a m + 0(a) . 
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Injecting this relation in the stationary energy balance for the zonal degrees of freedom (141) gives 

E z = 7zl x fdyE R (Yv&W) u\ + a z = a z + a m + 0(a) = 1 + 0(a). 



The barotropic equations (|5} are in units so that E = 1 in a stationary state: the above relation expresses the fact that, in the limit 
a <C 1, all the energy is concentrated in the zonal degrees of freedom: E z = E + 0(a). 



4.2.2 Enstrophy balance for the kinetic equation 

We conclude considering the enstrophy balance for the kinetic equation. As we will see below, the concentration property 
found for the energy does not hold for the enstrophy. The zonal and non-zonal enstrophy balances can be obtained with a very 
similar reasoning as the one done for the energy and, at leading order in a, one can use the full Fokker-Planck dl7l > or the 
approximated one P3l22t and obtain consistent results. Denoting Z z = E [A f@q%] the enstrophy of zonal degrees of freedom 
and aZ m = aE [A f@ afc] the non-zonal degrees of freedom one, we have 



= -2nLE 



= 2k1 x E 



J &yU"(y)(v [ i ) w m )(y) 



dyU"(y)(v ( i ] co„ 



(y) 



2Z z + 4tz 2 I x C z (0) : 



■2aZ m +47r z l x C m (0), 



(43) 



(44) 



where the first equation is the stationary enstrophy balance for Z z and the second one for Z m . The above equations refer explicitly 
to the case of the 2D Euler equations, but the generalization to the barotropic equations is straightforward. 

As will be discussed in next section, the enstrophy in the non-zonal degrees of freedom doesn't converge in the inertial limit: 
more precisely, aZ m = 0(1). Then, the enstrophy in the zonal degrees of freedom is 

Z, = - aZ m + 2k 2 l x (C z (0) + C m (0) )=I Z - aZ m ., 

with the total enstrophy injected by the forcing Iz- Then, by contrast with the zonal energy, the zonal enstrophy is not equal to 
the enstrophy injected plus correction of order a. There is no concentration of the enstrophy in the zonal degrees of freedom, 
and a non-vanishing part of the enstrophy injected is dissipated in the non-zonal fluctuations. 

Clearly there is no self-similar cascade in the problem considered here, energy goes directly from the fluctuations at any 
scales to the zonal flow through the effect of the advection by the zonal flow. These results are however in agreement with the 
phenomenology of a transfer of the energy to the largest scales while the excess enstrophy transferred to the smallest scales, 
however with dynamical processes that are non-local in Fourier space. This non-equilibrium transfer of energy to the largest 
scales is also consistent with predictions from equilibrium statistical mechanics which, roughly speaking, predicts that the most 
probable flow concentrates its energy at the largest possible scales. 



5 The Lyapunov equation in the inertial limit 

As discussed in section[3] it is essential to make sure that the Gaussian process corresponding to the inertial linearized evolution 
of non-zonal degrees of freedom close a base flow U, see Eq. d 1 3 b and l!24t . has a stationary distribution. We discuss this issue 
in this section. We consider the linear dynamics with stochastic forces Eq. l !22t . that we recall here for convenience 

^+Ll[a m ] = V2j lm , E[T) m (ri,f)T? m (r 2 ,/')] = C m (r l -r 2 )d(t-t'), (45) 

where 

L° W = U(y)^ + (-U"(y)+h'(y))^ (46) 
ox dx 

is the linearized equation close to the zonal flow U. 

Eq. < !45b describes a linear stochastic Gaussian process, or Ornstein-Uhlenbeck process. Thus, it is completely characterized 
by the two-points correlation function g(ri,Tz,t) = E [a> m (ri,r)»,„(r2,f)]. The evolution of g is given by the so-called Lyapunov 
equation, which is obtained by applying the Ito formula to ( 145 1 : 

§+L^+4 (2 V = 2C m . (47) 

We will prove that equation ( 147 ) has a asymptotic limit g°° for large time. This may seem paradoxical as we deal with a linearized 
dynamics with a stochastic force and no dissipation mechanism. We explain in this section that the Orr-mechanism (the effect 
of the shear through a non-normal linearized dynamics) acts as an effective dissipation. However this effects is not uniform 
on all observables. We will prove that g has a limit g°° in the sense of distribution, from which we will be able to prove that 
velocity-like observables have a limit. By contrast, observables involving the vorticity gradient will diverge. Moreover, if the 
kinetic energy contained in the non-zonal degrees of freedom converges, the enstrophy diverges. This non-uniformity for the 
convergence of the two-points correlations functions is also related to the fact that the convergences will be algebraic in time, 
rather than exponential. The statement that "the Gaussian process corresponding to the inertial linearized evolution close to a 
base flow U has a stationary distribution" must thus be understood with care: not all observable converge and the convergences 
are algebraic. 
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An observable like the Reynolds stress divergence involves both the velocity and the vorticity gradient. It is thus not obvious 
that it has an asymptotic limit. We will also prove that the long time limit of the Reynolds stress divergence — Ey ^(vm^ffl^)J 

and of its gradient — Ey |"^;(vm'a) m )j are actually well denned. The results in this section ensure that the asymptotic expansion 
performed in section|3]is well posed, at leading order. 

Because some quantities (such as the enstrophy) diverge in the inertial limit, it is of interest to understand how the Gaussian 
process is regularized by a small viscosity or linear friction. This corresponds to replace the operator Ly in Eq. fl47 b with the 
operator Ly defined in Eq. d!2b . Moreover, to be able to separate the effect of the viscosity and of the Rayleigh friction, we will 
introduce in the following the operators 

Ll [co m ] = L° [(Om) - vA co„, , (48) 

and 

L" [co m ] = L° [©„,] + aco,„ , (49) 

in which the superscript indicates which of the terms have been retained from the ordinal operator Ljj. 

For a = v = 0, the two-point correlation function g converges as a distribution. g(r,r',f) diverges point-wise only for value 
of y and y 1 such that U (y) = U(y'), for instance y = y'. We will prove that for small a this divergence is regularized on a universal 
way (independent on U) close to points (r,r') such that U(y) = U(y'), over a scale X a = » where k is the wavenumber 
of the meridional perturbation. We stress that only the local behavior close to the divergence point is universal. All quantities 
converge exponentially over a time scale 1 /a. 

We will argue that a small viscosity v also leads to a universal regularization of the divergence, over a scale A v = 

/ \l/3 

( ku\y) ) ' When both a and v are not equal to zero, the way the singularity is regularized depends on the ratio of the length 
scales X a and A v . More detailed results are discussed in this sections. 

Our analysis is based on the evaluation of the long time asymptotics for the deterministic inertial linearized dynamics (7). 
As far as we know, such theoretical results are currently available only for the linearized Euler equation [7 1, h = jS = 0. We will 
discuss the Lyapunov equation in the case jS ^ in a forthcoming paper. 



5.1 General discussion 

5.1.1 Solution of the Lyapunov equation from the solution of the deterministic linear equation 

In this first section, we show how to obtain a formal solution of the Lyapunov equation fl47b using the solution of the deterministic 
dynamics 4z-\-Ly with appropriate initial conditions. This discussion can be trivially extended to the cases in which the operator 
appearing in Eq. d45l > is L\j, L" or Lu, by simply replacing the operator L° with the appropriate one wherever it appears. 
We expand the force correlation function C,„ in Fourier series, we have 

C m (x,y) = £ c M cos(kx + ly). 

k>QJ 

We note that because C m is a correlation, it is a positive definite function. This explains why sin contribution are zero in this 
expansion. Moreover for all IgN* and I 6 Z, we have c>; > 0. The expression Q/cos(fct + ly) +c^_;COs(fcc— ly) is the most 
general positive definite function involving the Fourier components e lkx and e'*' or their complex conjugates. Here we have 
assumed that the correlation function is homogeneous (it depends only on x\ — xi and on y\ — yq). Its generalization to the case 
of an inhomogeneous force, for instance for the case of a channel would be straightforward. 

Because the Lyapunov equation is linear, the contribution of the effect of all forcing terms just add up 

g= 52 c ki8ki, (50) 

k>Q.l 

where g& is the solution of the Lyapunov equation l !47t with left hand side 2cos(foc + /y). By direct computation it is easy to 
check that ^ 

gki{ri,n,t) = I e-" L °f [e w ](x 1 ,y 1 )e-" L ?'[e^](x2,y2)dfi+C.C.. (51) 
Jo 

with &kl( x >y) = e'( kx+ly \ and where C.C. stands for the complex conjugate^. We note that e - ' 1 ^ [e#] is the solution at time fi of 
the deterministic linear dynamics d t +Ly with initial condition e^. 

Let us observe that from the solution of the Lyapunov equation, we can also easily obtain the evolution of E [<» m (ri,f) {Sco m ) (r2, 
where S is a linear operator. We have 

E[co,„(r u t) {Sco m )(r 2 ,t)}=S^[g} 

where is the linear operator S acting on the second variable. The typical operators S that we will use in the following are 
those necessary to obtain the stream function and the velocity from the vorticity. 

2 A similar formula can easily been deduced for any stochastic force of the form C(ri,T2) = f(Ti )/(i"2) from the explicit solution of the Gaussian 
process from the stochastic integral J ' e~'' L f [/] dW ti . Here the key point is that the Fourier basis diagonalize any translationally invariant correlation 
function. 
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5.1.2 Two explicitly solvable examples 



Before discussing the long-time behavior of the Gaussian process in Eq. d45b for a general zonal flow U, we discuss here two 
simple examples for which the deterministic linear dynamics d t + L^j can be solved analytically. The two examples are the 
perfect and the viscous advection by a linear shear. These two simple examples will put in evidence the mechanisms that will 
ensure the convergence of the long-time limit of two-points correlations in the inertial limit. 

For a linear shear this mechanism is the Orr mechanism: the transport of the vorticity along each streamline leads to a phase 
mixing in the computation of all integrated quantities, for instance the velocity. Then the velocity or the stream function decay 
algebraically for large times, the exponent for this decay being related to the singularity of the Laplacian green function. For 
more general profile U for which t/"(y) ^ 0, this shear mechanism still exist, but is also accompanied by global effects due to 
the fact that vorticity affect the velocity field globally. In many cases, those global effects are the dominant one |7 |. They will be 
taken into account in section [5721 



Perfect advection by a linear shear To treat an example that can be worked out explicitly, we consider the perturbation by a 
stochastic force of a linear shear U (y) = sy, which corresponds to set in Eq. ( 145 \ L° = s y^x- For sake of simplicity, we only treat 
here the case s > 0, the corresponding generalization to s < being trivial. The following discussion applies to flows that are 
periodic in the longitudinal direction x, with period 2kI x , either to the case of the domain & = [0,2kI x ) x (— °°,°°), or to flows 
in a zonal channel with walls at y = ±L. 

The Lyapunov equation we have to consider for this problem is 



■+L 0(1, e + L 0(2) e-2C 



for the vorticity- vorticity correlation function g(ri,r2,f) = E[<» m (ri,f)a> m (r2,?)]. For sake of simplicity, we consider the case 
where the forcing acts on a single wave vector 

Cm = £ J£±3. \j«M-*i) + il(n-yz) +cc ] = £ ( fc2 + /2 ) C os [k( Xl -x 2 )+l( yi - yi )] , (52) 

where C.C. stands for the complex conjugation of the first term. As explained in the previous section, contributions to the 
Lyapunov equation from other forcing modes just add up (see equation (150b). £ is the average energy input rate per unit of mass 
(unit m 2 s~ 3 ) (in this section and the following ones dealing with the case of a linear shear, co m is the non-zonal vorticity and has 
dimensions s , whereas in section dealing with the kinetic theory yfa(D m is the non-zonal vorticity). 
The deterministic evolution of the linearized dynamics with initial condition e'^ v+ '- v ' is 

e -'i" r e i(fct+fy)l _ e ik.x+ily-iksyr ^) 

as be easily checked. From i5 It . we have 

e (iA + l 2 ) e ik ( x t- x 2)+u(yi-y2) r ., , , -| 

g(ri,i*,0= V .., r- e-^-»)'-l +C.C, (54) 

4s -ik(yi-y 2 ) L J 

e(k 2 +l 2 ) 

foryi ^y 2 and g(xi ,y,x 2 ,y) = y 2 ' cos (k{x\ -x 2 ))t. 

This result readily shows that the square of the perturbation vorticity g(r, r) diverges proportionally to time t . This is expected 
as the average enstrophy input rate per unit of area is e (k 2 + 1 2 ) , and there is no dissipation mechanism. However, we also remark 
that for yi ^ y 2 the autocorrelation function is a fast oscillating function. As a consequence g will have a well defined limit in 
the sense of distributions: 

g"(ri,r 2 ) = lurig(ri,r2,0 (55) 



e r + l )_ Jk(xi-x 2 )+U(yi-y2) I 
As k 



x8{ yi -y 2 )-iPV 



yi -n 



+ C.C. 



where PV stands for the Cauchy Principal Value (we have used Plemelj formula, equation d94l > on page 129b). Equivalently, we 
have 



*"(ri,r 2 ) 



k 2 + l 



2s 



cosk(x\-x 2 ) sin[k(x\-x 2 ) + l(y\-y 2 )] ( 1 



n5 (y, -y 2 ) + 2±1 • i±iL FV 



yi -yi 



(56) 



The fact that the stationary vorticity-vorticity correlation function has a limit in the sense of distributions is a very important 
result. It means that every observable that can be obtained by integration of a smooth function over g°° will have a well defined 
stationary limit. Actually the formula above will be valid when integrated over any continuos function. For instance, we can 
use it to compute the velocity-vorticity or velocity-velocity correlation functions. Then all these quantity will have a definite 
stationary value. This is a remarkable fact, as we force continuously the perfect flow and no dissipation is present. Looking at 

the prefactor ^ - — '-, we remark that it is an injection rate e (k 2 + Z 2 ) (the injection rate of enstrophy per unit of area) divided 
by twice the shear. By analogy with the equivalent formula in classical Ornstein-Uhlenbeck process, we see that the shear s 
acts as an effective damping mechanism. The effect of the shear, called Orr mechanism, leads to phase mixing which acts as an 
effective dissipation for the physical quantities dominated by the large scales. All integrated quantities, for instance the kinetic 

energy, are proportional to v 2s — and are independent on the linear friction or viscosity at leading order. We also note that in 
the statistically stationary state the enstrophy is infinite. 
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As an example, we compute the vorticity-stream function stationary correlation function 

r(ri,r 2 ) = limE[(B(r 1 ,f)^(r2,0] (57) 

From Eq. d56t we obtain 

~cosk(x\-x 2 ) . f sin[k(xi -x 2 )+ 1 (y i - y 3 )] H k (y 2 -y 3 ) 



h (ri,r 2 ) =■ 



2s 



[ n H k{yi -y 2 ) + Pv j sm^i-^) + /( y i-y3)j^(y2-y3) dy3 

J k V1-V3 



(58) 



where Hj is the Green function for the Laplacian in the k! h sector, with the appropriate boundary conditions. From this last 
expression, and using v« = —dy/dy and vW = dy/dx it is easy to obtain the correlation function between the vorticity and 
the velocity field. We note that the Green function Hk{y\ —y 2 ) is not a smooth function: its derivative has a jump in y\ =y 2 . 
This implies that the correlation between the vorticity and the x component of the velocity is defined point- wise only foryi 7^ y 2 
or defined globally as a distribution. With analogous computations one can deduce other stationary two-points correlations: for 
example, the Reynolds stress divergence converges to a well defined function. 

In this section, we have discussed how the vorticity-vorticity correlation function has a stationary limit in the sense of 
distribution. When looking at this solution point-wise, it is singular for y\ = y 2 . In the following section, we explain how this 
singularity is regularized either by a linear friction or by viscosity. The discussion in this section was relying on the analytic 
solution of the linear dynamics close to a linear shear (Eq. (153l >). Such an analytical solution is not known for generic base 
flows, so that we need more refined techniques, as will be explained in section 15.21 We will obtain the same conclusion: the 
vorticity-vorticity correlation function converges as a distribution, and is regularized in a universal way. 

Advection by a linear shear: regularization by a linear friction We consider the same problem as in the last paragraph, but 
adding a linear friction. We will see how the singularity close to yy = y 2 is regularized by a linear friction. We solve 

dga a(l) T a(2) _ „ 

<~^u ga+^u ga — &~mi 

with Ljj = syd x — a and C m given by equation d52b . It is straightforward to observe that the evolution under the linearized 
dynamics L" is given by 

t -tL% r e '(fa+'v)l _ e ikx+ily-iksyt & -at 



With very similar computations to those of the last section, we can obtain the stationary value of the vorticity-vorticity correlation 
function 

g«(n,r 2 ) = limg a (n,ri,0 = £ (** + ^ jKn-vhtbi-yt) X F (y x - yi) +C .C., (60) 

t^°° 4s K Ts 

where 

y — iA 

We can thus observe that the function Fx (y) is a regularization of the Plemelj formula (equation ( I94t on page!29ll) on the length 
scale A . Indeed we have 

where the real part of Fx is even while the imaginary part is odd. Moreover the real part of Fx is a regularization of 7cS(y), 

and the imaginary part of Fx is a regularization of —PV(l/y), 

lim 3 [Fx (y)] = - lim . y = -PV (- 

where 9^ and 3 stand, respectively, for the real and the imaginary parts. We note that for y 2> A, lim^_5.o+ 9? [Fx (y)] = 0, and 
lim A ^o+3[F A (y)] = 1/y. 
We thus finally obtain 

g~( ri ,r 2 ) ~ t £ ) |cosfe(xi -x 2 ) Igl [f^(yi -y 2 )] (62) 

-sin[A;(xi -x 2 )+l{yx -y 2 )] ~3 [^(^i - 3^)] 

where one should observe that, because $H[F a (fcy)] decays sufficiently fast to zero for yS> a, the factor cos \k{x\ — x 2 ) +/(yi — y 2 ) 
has been replaced by cos [k{x\ —x 2 )]. Eq. l !62b is a regularization of the vorticity-vorticity correlation function found in (I56t , by 
the effect of a small linear friction a > 0. 

Quantities which were found to be divergent in the last section are now regularized by the presence of a small Rayleigh 
friction. For example the point-wise rms. non zonal enstrophy density 

% H] = gZ (0, 0) = £( ^ /2 ) . (63) 
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Advection by a linear shear in a viscous fluid In this paragraph, we study the regularization of the solution to the Lyapunov 
equation by a small viscosity. To treat an example that can be worked out explicitly, we consider the perturbation by a stochastic 
force of a linear shear U(y) = sy in a viscous flow. We consider the domain & = [0,2tzI x ) x (— °°,°°) and periodic boundary 
conditions in the x direction. 

We solve the Lyapunov equation 



2C»j 



for the vorticity-vorticity correlation function g(ri,r2,f) = ~R[(Q m (r\,t)(O m (r 2 ,t)\, where = syd x — vA,C„, is given by equation 
i and e the injection rate of energy per unit of mass (e has the dimensions m 2 s~ 3 ). 

The deterministic evolution of the linear dynamics 4j +Ly = with the initial condition e '( kx +b') j s given by 



,i(kx+ly] 



& ikx & i(l-skt)y - 



(64) 



This solution, first derived by Lord Kelvin, can be obtained by the method of characteristics. Alternatively one can directly 
check a-posteriori that (1641 is a solution to the deterministic equation. We can see that the second exponential in d64l > just gives 
the inertial evolution of the perturbation and the third one describes the effect of the viscosity. We thus see that the vorticity is 
damped to zero on the time-scale (l/i 2 fc 2 v)'/ 3 . This time scale is the typical time for an initial perturbation with longitudinal 
scale 1 jk to be stretched by the shear until it reaches a scale where it is damped by viscosity. 
We can now calculate the asymptotic solution to the Lyapunov equation. It is given by 



gv( r l> r 2) 



e{k 2 + l 2 



e -' L u 



v(l) 



-ikx2~—ily: 



e(k z + l 2 ) . k(xj _ X2)+U{yi _ yi} 1 



4s 



H, 



2v\l/3 



(yi-y 2 ) + C.C. 



where we have used Eq. d64t and the change of variable t\ = (2vk s ) /J f to write the second equality, and where 

A Jo 



(65) 



(66) 



In appendix|B]we prove that Hx (y) = Gx (y) + with 



Gx(y) 



>1 1 1 (3 
g-t'l-j't 



(67) 



We note that Hi(7) = GxiiAY) /n is one of the two Scorer's functions, that solve the differential equation d 2 Hi/dY 2 — FHi 
1 /n, and is related to the family of Airy functions. We stress that the asymptotic behavior for large y of the real parts of Gx and 
H x are different Si [Hx (y)] 



(k 2 + l 2 )^ whereas 9? [G A (y)] 



~ 2^3- but in any cases those are subdominant for small 



A (please see appendix [B] for more details). 

As explained in appendix |Bj Gx is a regularization of the distributions in Plemelj formula at the length scale X and has 
the dimension of the inverse of a length. The integral over y of the real part of Gx is K, consistently with the fact that it is a 
regularization of k8 (y) ; the imaginary part is a regularization of Cauchy Principal value of 1 /y, with 



We thus conclude that 



3[G A (y)] 



gv( r U*2) 



e (k 2 + l 2 ) f cos[fc(xi -x 2 ) 



(^)' /3 '«i 

sin [k(xi-x 2 )+l(yi-y 2 )} ( 



2s 



9? 



{if) 



(68) 



(69) 



G (|i)'/3(>'i-y2) 



Observe that, because 9? 



G (r)' /3 ^ 



decays sufficiently fast to zero for y S> 



[j^) 1 ^ 3 , the factor cos \k{x\ -x 2 )+l{y\ -y 2 )} 

^2v^l/3 



has been replaced by its value for y i = y 2 . The first term of Eq. d69t is a local contribution, for values of y i — y 2 of order ( ) 
whereas the second term is a global contribution. Whereas the local contribution is independent on / and depends on yi — y 2 
through the shape function G, the global contribution has a phase dependance through /(yi — y 2 ). 



Looking at the prefactor 



e(/t 2 +/ 2 



we remark that it is the ratio of the power input e (£: 2 + / 2 ) divided by the shear. By 



analogy with the equivalent formula in classical Ornstein-Uhlenbeck process, we see that the shear s acts as an effective damping 
mechanism. Integrated quantities are independent on the viscosity at leading order, for instance the kinetic energy will be 

- x 2l y\—y 2 ) is, independent on the viscosity at 



e(k 2 +l 2 ) 

proportional to v 2; — '-. Also the two-point vorticity correlation function gy(x\ 



leading order for values of (yi — y 2 ) much larger than 



'2v\l/3 
v ks 1 
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1 /3 

By contrast, the point wise value of the two-point vorticity correlation function for (yi — yi) ~ (jr) diverges for large 
Reynolds number. For instance 

e(k 2 + l 2 ) i s ^ i/3 

and one can note that j^i is a Reynolds number based on the local shear and the scale of the non- zonal perturbation. We observe 
that the enstrophy density Ey [co 2 ] is regularized by viscosity but diverges for large Reynolds number to the power 1 /3. 

We conclude by stating that physical quantities involving higher order derivatives will also be regularized by the viscosity 
and diverge with the Reynolds number. For instance the palinstrophy density will diverge as 



e{k 2 + l 2 ) /ks\ 2/3 , , ^ : 



Et/[(Vo) m ) 2 ] ~ C P . . . 

(§£) 1/3 {<i 2s V 2v / \2vk 2 

where Cp is a non-dimensional constant. 

In this section, we have discussed the case of a force at a longitudinal scale 1 /k and transverse scale 1 //. The conclusion for 
any other forces can easily be obtained by superposition of the contributions from all scales. The general conclusion is that in the 
limit of large Reynolds number 3> 1, the two-point correlation function converges as a distribution, and converges point- 

,5 ,' 

velocity- velocity and velocity-vorticity correlation functions have a limit independent on the Reynolds number, and the kinetic 
energy density (m 2 .s~ 2 ) is roughly proportional to j- s where e is the average energy input rate. 

We have made explicit computations only in the case of a linear shear U(y) = sy. Explicit computations are not easily done 
in more complex situations in the presence of viscosity. However we expect that for generic U, the singularities of the stationary 
solutions to the Lyapunov equations are regularized in a universal way, as argued in section 1271 



wise for values of y\ — y2 much larger than ( I where l/k m is the maximal scale for the forcing. As a consequence, the 



Advection by a linear shear in a viscous fluid with linear friction When both linear friction and viscosity are present, the analysis 
above can be easily generalized. The way the two points correlation function is regularized depends on the relative value of the 
two length scales \ a = ^2 and A v = (tt) • When X a \ v , the regularization is of a friction type and formula (1621 will be 
correct. When X a <A V , the regularization is of a viscous type and formula ( 169) will be correct. 

We stress that, whatever the values of the length scales X a and Ay, the real part of the regularizing function always decays 
proportionally to 1 /y 4 for small enough values of y and proportionally to 1 jy 2 for large enough values of y. The location of the 
crossover between these two behaviors depends on the values of the length scales X a and A v . A careful discussion of this issue 
is addressed in appendix iBl 

Here, we only point out that three cases are possible: (?) when X a 2> A v the regularization is of friction type and the crossover 
happens for y ~ X a \ (ii) when A a <C Xy(k 2 + 1 2 ) the regularization is of viscous type and the crossover happens for y ~ 
\/2/(k 2 + l 2 ); (Hi) when \ a <C A v but X a 2> Xl(k 2 + 1 2 ) the singularity is also of viscous type, but in this case the crossover 
happens for y ~ \J 2Ay / \ a . The emergence of this last length scale \J 2Xy / X a is due to the different multiplicative factors of 
the 1 jy 2 decay due to friction and in the 1 /y 4 decay due to viscosity (please see appendix IBl. 



5.2 Stationary solutions to the Lyapunov equation in the inertial limit 

In the previous section, we have considered a base flow with constant shear, for which analytical solution to the Lyapunov 
equation can be computed, and provides a qualitative understanding of its solution. We have concluded that it has a stationary 
solution, in the sense of distributions, even without dissipation. This solution diverges locally for y\ = yi, related to an infinite 
enstrophy. We have explained how those local divergences are regularized by a small linear friction or viscosity. The aim of 
this section is to prove that the same conclusions remain valid for any generic shear flow U(y), which is assumed to be linearly 

stable. We also prove that the Reynolds stress divergence — [(vi^Gin)] ar >d its gradient — ^(vm' ) a» m )j have finite values. 
Convergence results for other two-point correlations are also discussed. 

At a rough qualitative level, the reason why this is valid is the same as the one discussed in the previous section: the effect 
of the shear (Orr- mechanism). However this explanation can not be considered as satisfactory. The constant shear flow case 
has a zero gradient of vorticity, that's why the equation are so simple and can be solved analytically. Whenever the vorticity 
gradient is non-zero, the hydrodynamic problem becomes drastically different, coupling all fluctuations globally due to the long 
range interactions involved through the computation of the velocity. Any explanation based on local shear only is then doubtful. 
Moreover, most of jets in geophysical situations have points with zero shear U'(y) = 0. This is also a necessity for jets in doubly 
periodic geometries. Then the local shear effect can not be advocated. 

As already suggested, explicit analytical results are hopeless in this case. In order to prove the result, we rely on two main 
ingredients. First we use the fact that the stationary solution of the Lyapunov equation can be computed from solutions of the 
deterministic linearized dynamics, as expressed by formulas d50t and ( 151b . Second we prove that these formulas have limits for 
large times based on results on the asymptotic behavior of the linearized 2D Euler equations, discussed in the work J7]. At a 
qualitative level, the results of this paper show that the flow can be divided into areas dominated by the shear for which the Orr 
mechanism is responsible for a phase mixing leading to an effective dissipation. In other flow regions, for instance close to jet 
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extrema, where no shear is present, a global mechanism called vorticity depletion at the stationary streamlines wipes out any 
fluctuations, extremely rapidly. 

We discuss in this section only the case with no beta effect jS = 0, the case with beta effect will be discussed in a forthcoming 
publication. 



5.2.1 The Orr mechanism and vorticity depletion at the stationary streamlines 



In this section we summarize existing results on the large time asymptotics of the linearized Euler equations [7 |. We consider 
the linear deterministic advection with no dissipation. Because the corresponding linear operator Ly is not normal, a set of 
eigenfunctions spanning the whole Hilbert space on which Ly acts does not necessarily exist. Because U is stable, has no 
eigenmodes corresponding to exponential growth. Moreover it is a very common situation that the Euler operator has no 
modes at all (neither neutral nor stable nor unstable). A simple example for which this can be easily checked is the case of the 
linear shear treated in Section 15.1.21 We assume in the following that has no mode and /3 = 0. 

While the vorticity shows filaments at finer and finer scales when time increases, non-local averages of the vorticity (such 
as the stream-function or the velocity) converge to zero in the long-time limit. This relaxation mechanism with no dissipation is 
very general for advection equations and it has an analogous in plasma physics in the context of the Vlasov equation, where it 
is called Landau damping [53 1. This phenomenon was first studied in the context of the Euler equations in the case of a linear 
profile U(y) = sy in [ 14|. The work [7| was the first to consider the case in which the profile U(y) has stationary points y c such 
that U'(y c ) = 0, which is the generic case. 

We consider the deterministic linear dynamics d, a> + CO = with initial condition e lkx f(y). The solution is of the form 
G)(x,y,t) = e' fa <%(y,/). From Q, we know that0 

at(y,t) t ^e?(y)<r iku ® t . (70) 

We thus see that the vorticity oscillates on a finer and finer scale as the time goes on. By contrast to the behavior of the vorticity, 
any integral of the vorticity with a differentiable kernel decays to zero. For instance, the results for the x and y components of 
the velocity and for the stream function are: 

v * ^h^mm^t— ' (71) 

Vk [y,t) t^ik(U'(y)) 2 t 2 ' { ' 



Uy,t) -. 7^1- , 2 ■ (73) 



and 

°(y) Q- ikvl y)' 

In all the above formulas, higher order corrections are present and decay with higher powers in l/t. One should also observe 
that CO~(y) depends on the initial condition f(y). The asymptotic profile C0™(y) could be computed numerically, for instance 
from the resolvent of the operator An essential point is that C0~(y) has in general no local approximation, it is not a simple 
function of the local shear but depends on the whole profile U. 

These results have been proven for every shear flow U, also in the presence of stationary points y c such that U'(y c ) = 0. 
Moreover, it has been proved in [7 1 that at the stationary points co£°(y c ) = . This phenomenon has been called vorticity depletion 
at the stationary streamlines. It has been observed numerically that the extend of the area for which G)J°(y c ) ~ can be very 
large, up to half of the total domain, meaning that in a large part of the domain, the shear is not the explanation for the asymptotic 
decay. The formula for the vorticity Eq. (170) is valid for any y. The formulas for the velocity and stream functions are valid for 
any y ^ y c . Exactly at the specific point y = y c , the damping is still algebraic with preliminary explanation given in (7), but a 
complete theoretical prediction is not yet available. 

Using these results, we will study the stationary solutions of the Lyapunov equation with no dissipation. 



5.2.2 The vorticity auto-correlation function converges to a distribution in the perfect flow limit 

We prove now that the Lyapunov equation has a stationary solution, in the sense of distributions. The hypothesis are the same as 
in previous section: U is linearly stable, has no mode, and jS = 0. 

In order to understand the long time behavior of the vorticity-vorticity correlation function, we consider Eq. (15 It . We 
consider Eq. ( 1701 ). where c%(y,f)e' fa is the solution of the deterministic equation with initial condition s +ily (we note that Cfy 
then depends on / through the initial condition e y ). One can check that the integrand of the r.h.s. of this equation behaves for 
long time as 

OkbuhMfati) ~ ar(y l )&r(y2)^ imyi) - u{y2)h - (74) 

For yi and yi such that U(y\) / U{yi) a computation analogous to the one in section 15.1.21 shows that gu, and hence g°° 
diverge proportionally to time t . This is related to an infinite value of enstrophy. 

In the following, we write formulas for the case when U is a monotonic profile. Then each frequency correspond to a single 
streamline. In the opposite case, two streamlines may have the same frequency, and resonances between streamlines should be 

3 We prefer to use in this section the notation with a tilde to denote the solutions of the deterministic linear dynamics d,+lSj. 
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considered. The formula would then be more intricate, but the result can be easily obtained from ( 170b and the conclusion that 
the limit exists is still trufl From ( 174b . we get 

gu(r U T 2 ,t)= [ t & k (y l ,hy kxi &k(y2,h)e- ikX2 dh+C.C. (75) 
Jo 

where g r is a function which remains continuous in the t — > °° limit. We thus obtain, for the stationary vorticity-vorticity corre- 
lation function 

and conclude that the Lyapunov equation has a stationary solution understood as a distribution. 

Another quantity of interest is the velocity auto-correlation function and the kinetic energy density. We note that the velocity 
auto-correlation function is a quadratic quantity, that can be obtained directly as a linear transform from the vorticity-vorticity 
correlation function g, or from the deterministic solution to the linearized operator. From ( 150b . we see that the contributions from 
each forcing modes add up. Then 

By [v m (ri,f)-v m (r 2 ,*)] = £ C w- E «( r i> r 2)> 

k>0,l 

with 

£«(ri,r 2 )=e i ^-^) f v k (y u u).v* k (y 2 ,u)du + C.C., (77) 
Jo 

where y k e lkx is the velocity of the deterministic solution to the linearized equations d t 6) + L^a = 0, with initial condition 
gihc+ily^ as m tne p rev i ous section. Alternatively, E k i = V^'V^g";, where V = —V (.)] x e z is the linear operator giving 
the velocity from the vorticity and \^\ resp. are the operator V acting on the first, or second variable respectively. From 
J72I71I I. it is clear that Ey and thus the velocity autocorrelation functions have finite values in the inertial limit, even if no 
dissipation is present in this limit. 

We note that OcEu [\ m (r\ , t) ■ v m (rj , t)] is the non-zonal kinetic energy density (the kinetic energy contained in the non-zonal 
degrees of freedom). We thus conclude that in the limit of very small a, the non-zonal kinetic energy is proportional to a, and 
that its value can be estimated form the Lyapunov equation with a = 0. Those results are extremely important, as they prove 
that the scaling we have adopted all along this work in order to make an asymptotic expansion is self-consistent. 

5.2.3 Convergence of the Reynolds stress divergence and of its gradient 

The stationary value of the Reynolds stress divergence — ^(vm'a) m )(y)J and its gradient — Ej/ \^{v>n®in)(y)\ are central 
objects in the kinetic theory. Indeed, they enter in the final kinetic equation when written, respectively, for the evolution of the 
zonal velocity 05b or for the evolution of the zonal vorticity d34b . We have thus to prove that they are finite. 

We note that the Reynolds stress is a quadratic quantity that can be obtained directly as a linear transform from the vorticity- 
vorticity correlation function g. From l !50b . we see that the contributions from each forcing modes add up. Then 



F(y) = Eu [(v^co m )(y)] = £ c kl F kl (y) 



k>0J 



with 



F u (y) = Km a> k (y,u)v { k y) *(y,u)d u + C.C., (78) 

where a> k e' kx and v^e'^ are the deterministic solutions to the linearized equations d t a) + L°jS> = with initial condition e' kx+dy , 
as in previous section. 

Proving that the long-time limit of the integral in Eq. d78b converges is straightforward. Using Eq. (170) and (172K one easily 
realizes that the integrands behaves as l/u 2 for large u and thus the integral converges. One conclude that the Reynolds stress 
divergence is finite. 

The computation needed to prove that also the gradient of the Reynolds stress divergence T.k>oj c kidF k i/dy exists, has 
to be done with more care. Using ( 178b , ( 170b and (172b . one may be led to the conclusion that a contribution proportional to 
l/u exists, that would lead to a logarithmic divergence once integrated over time. This is actually not the case because the 
computation shows that the leading order contributions to the Reynolds stress cancels out exactly, as can be checked easily. As 
a consequence, we conclude that also the divergence of the Reynolds stress is finite. 

In order to reach this conclusion, we have used ( 170b and ( 172b . which have been established for those values of y that are 
not stationary for the profile U(y). For a stationary point y = y c no theoretical results are available for the asymptotic behavior 
of the velocity field. However, based on numerical evidences reported in (7), we discuss in Appendix [C] that the Reynolds 
stress divergence is finite also for y = y c . This will be also checked by direct numerical computations of the Lyapunov equation 
discussed below. 



We note that the difficult theoretical result related to those resonances is to establish j701 
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5.2.4 The effect of the small Rayleigh friction and viscosity on stationary solutions of the Lyapunov equation 

We briefly comment on the effect of the small linear friction, or of a small viscosity. From the discussions in sections 15.2. H and 
15.2.31 we have proved that the kinetic energy density and Reynolds stress have limit values independent on a or the viscosity. 
By contrast, the vorticity-vorticity correlation function g°°(ri,r2) diverges point-wise for yi = y 2 (more generally for any two 
points for which the streamlines have the same frequency), and that it has a well defined limit as a distribution. This also implies 
that the enstrophy contained in non-zonal degrees of freedom is virtually infinite for a = 0. We address how this is regularized 
for small values of a, or due to a small viscosity. 

We have already discussed this issue for the particular case of a linear shear U(y) = sy, in section [5. 1.21 both with and 
without viscosity. For the general case of any stable base flow U and for a linear friction, we will conclude below that the same 
conclusion as for the linear shear holds. More precisely, we can prove that in the vicinity of y\ = yi, the singular behavior is 
regularized in a universal way (with shape functions independent on U) over a scale a/kU'(y). Moreover we can prove that 
E[/ [cfl^(r)] diverges proportionally to 1/of, such that the non-zonal enstrophy density CtEy [ffl^(r)] has a finite limit (this could 
have been expected in order to balance the finite enstrophy input rate provided by the stochastic force). 

When viscous dissipation is present, we think that it is still true that the divergence of g has a universal regularization, and 
we think that the scaling for the divergence and enstrophy found for the linear shear are also valid for a generic velocity profile 
U. We do not prove this result, but discuss why this is plausible at the end of this section. 

In order to prove those results, we need to consider the operator L" = L^ + a instead of the inertial for the analysis of the 
Lyapunov equation. The effect of replacing Ly by L" in section [5~. 1 . 1 [ corresponds simply to add the exponential damping e~ 20 "' 
to the integrands of expressions like the one in Eq. d51t . This observation gives us the possibility of understanding how diverging 
quantities are actually regularized by the presence of a small friction for any zonal flow U. Straightforward computation then 
leads to a regularization of the inertial result (1761 as 

gk ' {rur2) = \ik(U( yi )-U(y 2 )) + 2a +8a{yuy2 -' ^j 6 +CC " (?9) 

where, as before, g r a has a finite limit when a — > 0. In the small a limit and for yi ~ y2, we have 

|2 

)_ 

r<i.yi~^ _ l v " \kU'(yi 

1 



*2(ri,ii) „ ~ 2(cos*(x 1 -* 2 ) ^k' U '. 9i 



+j[cCyi)«r(>'2)e^-^ 



F 2a (yi -y 2 ) 
)| 



(80) 



H/'Cyi 



f 2a , (yi -yi) 



+ g« 2 (ri,r 2 ) 



where g'u has a finite limit when a — > 0. We also recall that if U'(y) = 0, then cB£°(y) = 0. Comparing d62t with l !80b , we 
conclude that the vorticity-vorticity correlation function is regularized in a universal way close to yi = y%. 

The fact that the Ef/ [co 2 (r)] diverges proportionally to 1 /a, is obtained by observing that Eq. d74t should be replaced by 

®k(yi,ti)c%(y2,h) ~^_>- C(yi)C*fc)e-*[ £/ ^'- f/ feM f '- 2a " . (8i) 

We thus have 

E£/ H(r)] = I c„g5(r,r) = +g r a (y,y,°°) (82) 



We have not worked out the case with a small viscosity. This is technically possible. This would require first to describe 
precisely the large time asymptotics of the linearized 2D Navier-Stokes equations. Such results should involve the classical 
literature on the study of viscous critical layers [22], though matched asymptotics expansions between what happens close to 

the critical layers and away from it. Those classical results exhibit naturally the scale and the family of Airy functions. 

It should be possible to extend those classical deterministic results to the Lyapunov equation, using the relation between the 
deterministic solutions and the solutions to the Lyapunov equation discussed at the beginning of this section. For this reason, we 
guess that the regularization obtained for the linear shear is universal and should explain the regularization for generic velocity 
profile U. One important technical drawback is that we guess that the classical results for the deterministic equation do not yet 
exist yet for the case where U has stationary streamlines (U'(y c ) = 0) and when several streamlines have the same frequency. 



5.3 Numerical solutions of the Lyapunov equation 

In the previous section we have proved that the Lyapunov equation for the linearized Euler dynamics has a stationary solution. 
In this section we compute numerically this stationary solution. 

For the Lyapunov equation with viscosity, or hyperviscosity, assuming a stable base flow, we expect the stationary solution 
to be well behaved. Then a natural way to solve the Lyapunov equation would be to discretize the linear operator L^ 1 ' and 
to directly solve the approximate dynamics. This is the traditional way and such a technique has for instance been used in most 
of previous works using SSST-CE2 [1 26 25. 50. 67 58 1. 

However, we are here specifically interested in the inertial limit. Then one could solve the Lyapunov equation for finite 
values of a and v and then study the asymptotic behavior of the results when these parameters go to zero. While feasible, this 
route seems extremely difficult, as the numerical discretization would have to be increased as v goes to zero. We have then 
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chosen to try another route in the following. We will make a direct numerical computation of F, the divergence of the Reynolds 
stress which appears in the kinetic equation d35t . which we have proved to be well behaved. In order to do that, we first establish 
an integral equation verified by h k i, the vorticity-stream function stationary correlation function, which will be precisely defined 
below. As we will show, this integral equation does not suffer from the instability in the small viscosity limit that would hinder 
a partial differential equation approach, and can be solved even if v = 0. 

Using this technique, we illustrate the theoretical results of last section: F has a finite limit for the zero friction limit a — > 0, 
the energy in the non-zonal degrees of freedom has a finite limit, and the energy dissipation rate for the non-zonal degrees of 
freedom goes to zero when a goes to zero. We also check the analyticity breaking, and its order, for the correlation function 
huiyuyi) at the points yi = y 2 . 

5.3.1 Computing stationary solution of the Lyapunov equation from an integral equation 

By definition, the stationary vorticity autocorrelation function gT, is the solution of the stationary Lyapunov equation (with linear 
friction and no viscosity) 

L" 8u+ L u 8kl = 2cos(£(xi -x 2 )+l{yi -yi))- 
Equivalently, we have g kl {x\,x 2 ,y\,y 2 ) = ^ Xl ~ X2 ^ gkliy 1J2) +C.C. with 

CW<!l£«=e<''(>'>-^ (83) 

where 

L® k = ikU - ik (U" -P)A- l +a. 

(2) (2) 

We define %; such that g k i = A k hkh where Ai is the Laplacian operator with respect to the second variable. The drift term in 
Eq. d33b can be directly computed from h k i, with 

F(y) = E u [(v^a m )(y)] = £ c kl F kl (y), 

where 



k>0,I 



F kl (y) = 2k3[h kl (y,y)}. 
In terms of h k \, the stationary Lyapunov equation d83t becomes 

Al), , s ik(U"( yi ) -P)h* kl (y2, yi ) -ik(U"(y 2 ) ~ P) h kl ( yi ,yi) W^'^ 

\ hl{yi,y2) = 77777 777 rr ; r~ . (84) 

ik(U(y\)-U(y 2 ))+2a 

We note that this equation is not a differential equation for hu(yi,y 2 ) for each y\ fixed, as it involves h* kl {y2,yi). 

(2) 

In the previous section, we have shown that g^, and thus g^i = A k h k i diverge point-wise for yi — > y2 in the limit a — > 0. 



This is related to the vanishing of the denominator for a = and yi = y2 in equation d84l >. On the other hand, it can be proved, 
with a very similar reasoning to the one used in section 15.2.3 1 for F k j, that h k \ is well-defined as a function, even in the limit 
a — > 0. Thus, we chose to turn Eq. d84b to an integral equation. Inverting the Laplacian operator using the Green function H k 
(see equation j58b ). we obtain 

hkl{y ^ = -k.l u M -u { y 2 )-^ (85) 

+ f ru f^f? «a [(U"(yi)-P)hh(y 2 ,y l )-(u"(y' 2 )-P)h kl (y u y' 2 )]dy 2 . 
J U(yi)-U(y 2 )- — 

The generalization to the 2D barotropic equation upon a topography is straightforward, replacing jS by h'(y) in d85 b - The 
main advantage of this equation is that it involves only well-behaved functions, even in the limit of no dissipation a — > 0. 
Indeed, in this limit, the integrals converge to their Cauchy principal values. Moreover, the fact that it doesn't involve any space 
derivative will make it easy to solve numerically, as discussed in next section. 



5.3.2 Numerical implementation 



In order to numerically compute solutions of Eq. (1851 . we chose an iterative scheme. We compute the sequence {h^} N> Q with 

hN+i =S + T[h N ], 

where S is the first term in the right hand side of Eq. (185) and T is the integral operator of the second term. If this sequence 
converges, then h k i = lim^^co^. 

We note that we have not been able to establish conditions for which T is contracting. As a consequence, the convergence 
of the algorithm is not guaranteed, and we will establish the convergence empirically on a case by case basis. More precisely, 
the convergence of the iterative algorithm is checked by plotting log 1 1 /zjv — ^iv-1 II as a function of the number of iterations 
where . is the L 2 norm. We have recently devised a different algorithm that systematically converges and allows to compute 
stationary solution of the Lyapunov equation for zero viscosity. This much more robust and versatile algorithm will be presented 
in a following paper. 



Kinetic theory of jet dynamics in the stochastic barotropic and 2D Navier-Stokes equations 



23 




The computation of integrals of the form / u^J^^f^ia/k ^ re 1 u i res a particular attention. Indeed, the singularity of the 
denominator at the points such that £/(/) = U(y) can be the source of important numerical errors, and we find that the result 
strongly depends on the resolution if it is not precise enough. The resolution required to get robust results is easily understood: 
in order to avoid numerical errors, the denominator must satisfy \U(y) — U(y±Ay)\ <S for a discretization step Ay. More 
precisely, for sufficiently small Ay, we have U(y) — U(y±Ay) ~ ±U'{y)Ay, so that the condition becomes Ay <§C jjj^y- The 
numerical results confirm this scaling. This is an important remark, because the iterative algorithm may converge and yet give a 
wrong result if the condition Ay <§; is not respected. 

J. 3.3 Reynolds stress and stationary correlation functions 

The properties of F and discussed at a theoretical level in sections |5TI and I5T21 are now illustrated in two examples. 

Parabolic profile We consider the case of the 2D Euler equations (/3 = h = 0) in the channel (x,y) e [0,2k) x [—1,1] with 
rigid boundary conditions \jf m (x,±l) =0, \ m (x + 27i,y) = \ m (x,y), and with a parabolic base flow U{y) = A(y + 2) 2 — Uq, 
where the constants A and Uq are chosen so that the total energy is 1 and the total momentum is 0. This flow has no inflection 
point, whence its linear stability by direct application of Rayleigh's inflection point theorem 1221 . Moreover, this flow has no 
stationary points (yo such that U'(yo) = 0). We will consider next the case of a cosine profile, which has stationary points. 
We moreover chose here, for sake of simplicity, to force only one mode k = I = 1 . This corresponds to the forcing correlation 
function C(r) = c\\ cos(jc + ;y), with en = 4.72 a normalisation coefficient. 

We have numerically computed the solution of d85b with the iterative scheme previously explained. As mentioned above, 
the necessary resolution to use depends on the value of a; in the present case, it ranges from Ay = 1 /60 for the largest value 
a = 0.1 to Ay = 1 /300 for the smallest one, a = 0.005. 

FigureQ]shows the Reynolds stress divergence F(y) = levied [hici(y,y)\ and it clearly illustrates the convergence of F for a — > 0, 
as theoretically described in section I5T21 We also note that the Reynolds stress divergence and the base flow profile U have the 
same sign, except in a small region near y = 0. Looking at the evolution equation for U ( 1371 ) 

d,U = aF-2aU + a£,, (86) 

this observation implies that the Reynolds stress is actually forcing the zonal flow. 

Others correlation functions can be very easily computed from and the the analyticity properties that we have predicted 
can be directly checked. For instance, the vorticity auto-correlation function g^i(xi,X2,yi,yz) = e'^* 1- * 2 -' /VP' hki(yi,y2) +C.C 
is represented in Figure [2] We clearly see the divergence at the point y\ = yi = when a — > 0. Moreover, we recover the 
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Fig. 2 Real and imaginary parts of the k-th Fourier component of the stationary vorticity auto-correlation function A k h^i(y,y2i\y 2 =o m me case of 
a parabolic base profile in a channel geometry, with k = I = 1 and different values of the friction coefficient a. The plots clearly show the expected 
divergence at y = 0. 




universal shape of A^'hf,i(yi,y2) near this divergence, as expected from equation ( 179) : for ^ -C y <C 1 , where sq = U'(0) is the 
local shear at y = 0, we have 



,(2),..,,, n ,l 2a 3 r (2) i ^ -fc(t/(v)-t/(0)) 

fc 2 (t/(y)-[/(0)) 2 +4a 2 ' L 4 -I F((/(y)-y(0)) 2 +4a 2 ' 



The comparison between this theoretical prediction and the numerical results is shown in Figure [3] 

Cosine profile The second example we consider is the zonal base flow U (y) = cosy in the domain (x,y) = [0, 2nl x ) x [—11, Jr) 
with periodic boundary conditions, which is usually referred to as the Kolmogorov flow [7 |. This flow is stable and the linearized 
operator associated to this flow has no normal modes for aspect ratio l x < 1 (7). We choose the parameters l x = 0.5, k = 2, / = 0, 
corresponding to the forcing correlation function C(r) = C2ocos(2x), with a normalisation coefficient cjo = 1.29. The Reynolds 
stress divergence F is plotted in figure|4] It converges to a smooth function in the inertial limit. For all points such that U'(y) ^ 0, 
this was expected from the theoretical results of section [5] We note that we have also a convergence of F(y) to a finite limit at 
the stationary points y = 0, K, as discussed at the end of section 15.2.3 1 We observe that the Reynolds stress is forcing the flow 
except in some regions around the zeros of U, like in the case of the parabolic zonal flow. 

5.3.4 Conservation laws 

In the channel geometry considered, the linear momentum is always zero, P x = J dyv x = 0. This implies the constraint / dy U(y) = 
0. At the level of the kinetic equation j86b . this implies that the integral of the divergence of the Reynolds stress is zero: 
/ dyj [hki(y,y)] = 0, which is also a trivial consequence of the definition of ha. This constraint is fulfilled by the numerical 
results in Figure [T] within an error of order 10 ~ 8 . 
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We now discuss the energy and enstrophy balance derived in sections|4]and|5] The theoretical results predict a vanishing of 
the non-zonal energy in the inertial limit, and a divergence of the non-zonal enstrophy as I /a. Figure [5] shows the three rates of 
energy involved in (142 ) as a function of I /a: the power injected by the forcing (by definition equal to 1) 

a m = -2n\(A- l C m )(0) = l, 

the power transferred to the zonal flow 

P, = 7il x fdyEu [(iti ] co m ) (y)] U(y) = nl x J Ay2kc kl 3[h kl {y,y)]U{y) , 

and the rate of energy dissipated in the non-zonal degrees of freedom 

P diss = aE m = -2a J Ay9K[h a (y,y)] . 

It illustrates the fundamental property that in the inertial limit, all the energy injected in the fluctuations is transferred to the zonal 
flow, as expected from the discussion in section|4] We note that the energy balance (142 ) is verified in our numerical solutions up 
to errors of order 1CP 4 . 

Enstrophy has a very different behavior. From the general results of [5] we know that the dissipated enstrophy is not expected 
to go to zero as a — > 0, so the transferred enstrophy is not expected to be close to the injection rate. In the case of a parabolic 
profile U(y) = A(y + 2) 2 — Uq, this property is even more dramatic as the transferred enstrophy vanishes, it is indeed 



Z, 



-MM 



dyU"(y) 



Vm (0„ 



(>') 



= -Kl x / dy2A • 2kc kl 3 [/%(y,y)] = 



which is zero because it is proportional to the conservation of the momentum. Whatever the value of a, no enstrophy is trans- 
ferred to the zonal flow and all the enstrophy injected by the stochastic forcing is dissipated in the fluctuations. When it comes 
to the numerical results, checking that the enstrophy balance is fulfilled amounts at checking the conservation of the momentum 
f dy 3 [hki(y,y)} = 0, which has already been discussed. 
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Fig. 5 Stationary energy balance for the fluctuations, in the case of a parabolic base profile in a channel geometry, with k = I = 1 and different values 
of the friction coefficient a. We see that in the inertial limit a -C 1, all the energy injected in the fluctuations is transferred to the zonal flow, while the 
energy dissipated in the fluctuations vanishes. 



6 Attractors of the slow dynamics and multistability of atmosphere jets 

In section[3j we have obtained that the jet has a slow evolution whose dynamics is described by equation d37t . reported here for 
convenience: 

^=F[U] + ^ (88) 

where is Gaussian process with autocorrelations given by (136) . In most known situations, there is no noise acting on the largest 
scales. Then the correlation function of <§, see Eq. J36t . is of an order a smaller than F [U]. In the limit we consider, we thus 
expect the noise to be small. Then one may wonder how important is its effect for practical purposes. We discuss this issue here, 
with emphasis on the case where the deterministic equation ^ = F [U] has more than one attractor. 

First it is useful to note that the deterministic equation ^ =F[U] is unable to describe the statistics of the small fluctuations 
close to a an attractor Uq. A first very interesting result that can be derived from (1881 is the statistics for the Gaussian fluctuations 
of the jet close to its most probable value. 

If we now assume that the deterministic equation 4^ = F [U] has more than one attractor, for instance Uq and U\ , the 
small noise | becomes essential in order to describe the relative probability of the two attractors and the transition probabilities 
between them. These can be computed numerically from d88t by means of large deviation theory or, equivalently using a path 
integral representations of the stochastic process (188) and instanton theory. One may wonder if such bistability or multistability 
situations actually exist. 

Multistability and phase transitions are actually very common for turbulent flows (for instance, paths of the Kuroshio current 
[64 1, atmospheric flows 1711 , Earth's magnetic field reversal and MHD experiments 0, two-dimensional turbulence simulations 
and experiments [66 8 45], and three-dimensional flows [61 1 show this kind of behavior). Figure [6] shows random transitions 
in the 2D Navier-Stokes equations HI. These transitions may in principle be described in the kinetic approach developed in the 
present paper; however, this would require an extension of the work to cases where non zonal attractors are also considered. This 
is in principle possible, but the Lyapunov equation is then more tricky. We postpone such an issue for future works. 

Remaining in the zonal framework assumed in this paper, we refer to works by N. Bakas, N. Constantinou, B. Farrell, and 
P. Ioannou 1 26 , 25 , 19] for existence of bistability. In those three works the possibility for multiple attractors for SSST dynamics 
(the deterministic part of Eq. d88t ) is discussed. Figure [7] shows such a case, where the evolution of SSST dynamics has been 
computed for some fixed parameters and fixed force spectrum. The two spatiotemporal diagrams of the zonally averaged velocity 
profile U obtained from the barotropic equations, and the two asymptotic profiles U obtained both from SSST dynamics and the 
barotropic equations show that for this set of parameters the deterministic part of the dynamics has at least two attractors, and 
that these attractors are observed in the barotropic equations. In order to compute the relative probability of the attractors and 
the transition probabilities from one to the other, one can rely on the Fokker-Planck equation d33l > or on the stochastic dynamics 
, These transitions involve states that are arbitrarily far from the attractors Uq and U\ . 
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Fig. 6 Figure taken from [8 ] showing rare transitions (illustrated by the Fourier component of the largest y mode) between two large scale attractors of 
the periodic 2D Navier-Stokes equations. The system spends the majority of its time close to the vortex dipole and zonal flow configurations. 
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Fig. 7 Evolution of the zonally averaged velocity profile as obtained from non-linear simulations (left) and comparison of the stationary velocity profile 
obtained from non-linear simulations and from the deterministic part of the kinetic equation (88), or equivalently SSST [25. 19 26|. The upper and 
lower pictures are obtained for the same values of the physical parameters but with different initial conditions. The figure shows that for a given set of 
parameters it can converge towards two attractors, Uq and U\. Courtesy N. Constantinou. 



7 Conclusion 



In this paper, we have considered the kinetic theory of the equations for a barotropic flow with beta effect or topography and with 
stochastic forces, in the limit of weak forces and dissipation a — > 0. We have formally computed the leading order dynamics, 
using the framework of stochastic averaging for this adiabatic process. At leading order, we have obtained a Fokker-Planck 
equation, or equivalently a stochastic dynamics, that describes the slow jet dynamics d37t . We have discussed that, at leading 
order, the deterministic part of this dynamics are the equations known as Stochastic Structural Stability Theory or, equivalently 
what can be obtained from a second order closure theory. This provides a strong support for those equations if one wants 
to compute the attractors for the jet dynamics, when a is sufficiently small. This result also shows that those models do not 
describe completely Gaussian fluctuations close to these attractors. Moreover, in phenomena for which large deviations beyond 
Gaussian become essential, a kinetic approach is still relevant, and the full kinetic equation l !37t should be considered. 
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We have also proven that at leading order, quasilinear dynamics will lead to the same kinetic equation 07t in the limit of 
weak forces and dissipation. We stress however that this will be wrong at next order. 

Our main hypothesis is that the zonal jets with profile U (y) are stable and have no non-zonal neutral modes. Whereas the 
stability is a crucial hypothesis, the non-neutral mode hypothesis is in principle not essential. However, if this hypothesis would 
not be fulfilled, a different kinetic equation would be obtained. 

For the kinetic theory to be meaningful, the base flows have to be attractors of the inertial dynamics. This is actually the 
case for some of zonal jets. However non-zonal attractors also exist, for instance dipoles for the 2D Euler equations [ 18 8 |. The 
kinetic theory could be extended in principle to non-zonal attractors. It would be easy to write phenomenologically the kinetic 
equation, just as was done for zonal jets. We stress however that proving the self-consistency of the theory, which is the key 
point from a theoretical point of view, would be much more tricky. This requires the analysis of the asymptotic solutions of the 
Lyapunov equation, which would be technically more difficult for non-zonal attractors than for zonal jets. 

A number of interesting issues have not been addressed in this work and should be considered in the future. Whereas the 
limit a — > is the appropriate one in order to set up a clear theoretical framework, an essentially question arises from a practical 
point of view: how small a should be? How does the answer depend on the force spectrum? These issues could certainly be 
addressed both numerically and theoretically. 

We believe that our theoretical approach is the correct way to describe the statistics of the largest scale of the flow. However, 
for small but finite a, especially if the forcing has a narrow band spectrum, it is likely that non-linear effects would become 
important. In some regimes, we expect to see together with a quasi-linear dynamics at the largest scales, either enstrophy 
cascade, or inverse energy cascade, or both, at scales much smaller than a typical jet scale. Those are described by Kraichnan 
type cascades for the 2D-Navier-Stokes equations [4], or by zonostrophic turbulence if jS / I31II30I . In such a case, how 
would our large scale theory relate with smaller scale cascade pictures for finite a? What would be the crossover scale and how 
should it depend on a? While answer using dimensional analysis can be given easily, a more detailed theoretical study would 
be very interesting. Another related issue is the effects of the small scale part of the force spectrum. Allover the paper, we have 
assumed that the spatial correlation C is smooth and that its spectrum converges sufficiently fast for large wave numbers k in 
order to insure the convergence of the the sum of the contributions of all sector k to the Lyapunov equation. Could the required 
hypothesis be made more explicit? 

We have made explicit computations up to leading order only. It would however be extremely interesting to compute next 
order corrections. This may be important in order to have very precise predictions for small but finite a and also to give a clear 
indication of the order of magnitude of a required for the theory to be trusted. Undoubtedly the answer will depend much on 
the force spectrum. We stress that the next order corrections will be different from an approach using a cumulant expansion and 
that it should not be hindered with the usual inconsistency observed in cumulant expansions. 

Another important issue is the validity of our approach for more complex dynamics, for instance for the class of quasi- 
geostrophic dynamics. At a formal level, the same approach can be applied. However, the dynamics is expected to be quite 
different. For instance, whereas the fact that there is no non-zonal neutral modes is probably generic for barotropic flows, it is 
probably exceptional for equivalent barotropic quasi-geostrophic dynamics as soon as the Rossby deformation radius is smaller 
than the domain size. Also if one consider a deterministic force in addition to the stochastic one, which is essential for instance 
for Earth jet dynamics, then the same formalism can be applied, but the phenomenology is expected to be quite different. This 
calls for new specific studies. 

Finally, could the validity of the kinetic theory be proven in a mathematical sense? Our results on the Lyapunov equation 
can be considered as mathematically rigorous. Their consequences for the self-consistency of the asymptotic approach too, 
for instance that the Reynolds stress and other physical quantities computed from the Ornstein-Uhlenbeck process are finite. 
However, in order to claim for a full mathematical result for the kinetic approach, one would need to control more precisely the 
asymptotic expansion. 
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A Explicit computation of the operators in the zonal Fokker-Planck equation 

We compute here each term of the equation for the slowly varying part of the PDF B2t . which we report: 

^ = \a»X z + a 3 ' 2 ,i?3? z jT dt's/^jfn + a 2 »^ z ^°°d/'e'' r " [(1 - ,9>)^ z + t'^ 2 ] \p s + ff (a 5/2 ) . 

This will lead to the final reduced Fokker-Planck equation for the zonal flow PDF R, {33). 
- The first term in equation j891 gives 



}[q z ,(O m ] j 



Ay l 



dy 



(v^W)] iy\) + (O z {yi)-^Aa z {yi)\R+ J AyzC z iyi,y 2 ) 



5R 



The next term is exactly zero: we have 

X,P S = J dr, 



5<D,„(ri) 



[{v m .Va) m (T l ) - {\ m .Va,„(r 1 )))G[q : ,(O m ]R[q z }} . 



(89) 



(90) 
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We remark that the action of the functional derivative Sol ° will produce factors which are either linear or cubic in the variable o,„. When 
applying S? z , those terms will be multiplied by constant or quadratic factors. The projection will then lead to the computation of odd moments 
of the centered Gaussian distribution ((25), which are all null. As a consequence, we conclude that 

' z e TS "°J%P, = 0. 



(91) 

This is a very important result: at first-order, the corrections to the quasi-linear dynamics due to the full nonlinearity exactly vanish. This is a 
partial explanation to the success of the quasi-linear dynamics: the non-linear interactions between the non-zonal degrees of freedom (eddy-eddy 
interaction) don't have any influence on the statistics of the zonal flow until the order a 2 . 
Terms of order 0(a 2 ) in equation {89) give non-trivial contributions. First, we have 



; j dridr2 
5 



(v„,.VQ),„(r 2 )-(v m .V(B m (r2)))... 



Sra m (r 2 ) 

(v,„. Vq)„, (n ) - (v,„.Vo,„ (n)» G [q z , raw] r [&]] \ . 



leading to the computation of even moments of the Gaussian process, so the preceding remark don't apply here. Writing 

S%P, = M[q z ,W m }G[q z ,m m }R[q z ], 

we get 



+ 



G[q z ,a>,„] j dvi 
G[q z ,m m ] J dyi 



5 



Et/ {K> m m )( yi , t)M[q z , m m ] (0)]R[q z 
v 



o ; (yi) - -A(O z ( yi )j E v [M[q z ,(O m ]]R[q z ]... 



Sq z (yi 

+ 1 dnCt(y u n)gA^ (% [M[q z ,(0,„]]R[q z ]) 

Like for equation i29) , we have [M] = 0, so the two last terms in the above expression vanish. Using Liouville's theorem 

5 



dr- 



■ (v,„.Vco m (r)) = 



Sco,„(r) 

and the relation ( (/ — (/) ) g) = (f (g — (g))) , we find the expression of M 

M[q z ,a>„,]=- I dr!dr 2 (v„,.Vra,„(r 2 ) - (v m .Vco„,(r 2 ))) (v,„.V(B„,(ri ) - (v m .Ve&.(n )))g~' (ri ,r 2 ) 

r r I 2 

+ / drdr'(v m .V<B„,(r)-(v m .V(B m (r)))g- 1 (r,r')ffl„,(r') 

+ ydr,dr 2 dr (v m .Va),„(r 2 )-(v,„.VcB m (r 2 )))G(r 1 ,r 2 ).V ri (ffl,„(ri)) ... 

x ^-g- l {r l ,r)mjr) + (g- l (r l ,r)) ri co m (r)] 
+ / drdr 2 (v m .Vco m (r 2 ) - (v„,.Va>,„(r 2 ))) ... 



x -v m (r 2 ).V r ,g I (r 2 ,r)fij,„(r)+vS' ) (r 2 )a 3 , 2 (g 1 (r 2 , r) ) w„, (r) 



withg = g°°. 

In the last term, we have 



(l-^)J? z P s = G[q z ,w m ]J 



dv. 



Sq z {yi) [dy 



thus 



f d) ' 2 8q S (yi) { ( l!y~2 ( ^ C °'" ) ^ 2 ^ + ^ ) ~ ^ 4 ^ ^ ^ + f dy3 C ' & ' y3 ) 



Sq : b'3 



}[q z ,a>„] j 



dy, ■ 



d 
d7. 



Using the definition of the correlation and covariance | |261 and (27) we conclude that 
G[q z ,a> m ] / dy 2 dyi 



v^ffiw)-^ [(v2 (B,„)]) (yi)R j. 



~\2 

{Ec [[(v^ccvn) (yi,T) (vWffl m ) (y 2 ,0)]] Rfe]} . 



(92) 



(93) 



B Linear friction and viscous regularization of Sokhotskyi-Plemelj formula 

The Sokhotskyi-Plemelj formula adapted to our notations is 

/ e-" y dt = lim - izS (y) -iPv( - ) . (94) 

Jo a^a+ y-ia \y J 

In this Appendix, we consider different regularizations of the above equation that correspond to a viscous regularization or to the combination of the 
regularizations given by viscosity and by linear friction . We thus deduce few useful technical results. 
In the following, we have to deal with the large y behavior of oscillatory integrals of the form 



30 



Freddy Bouchet et al. 



r e -yf(t)dt. 

Jo 



In this respect, we will use the well known result that 



and 



/ e- yt f(t)dt =-Y V 



H f {2n) {{)) 
'2/1, ,2n+l 



(95) 



(96) 



(97) 



//=0 ^ 

The above formula holds for any function / which is (2N + 1) times differentiable, such that lim [/'"'(')] = f° r an Y n ^ anQ sucn that 
Jq dte' y, f^ (t) is finite for any n < 2N + 1 . They can be found by successive part integrations. 

Viscous regularization. We consider the regularization of Plemelj formula that correspond to the effect of viscosity. We consider 

H X M = I f° e -'x«- x2 (* 2 +' 2 )»+^"M» 3 d». (98) 
A ,/o 

We note that the real part of Hi is an even function and the imaginary part of Hi is an odd function. Moreover, for every fixed F, we have Hi (XY) 

• ■ 

G x {XY) + &{X), where G x is defined in Eq. <67l 

We can also compute the integral of Hi using the fact that 



fa (y)dy = ±- Jl dy /" e A>-^< 2 )<^> 2 -\> 3 dt , 



(99) 



which can be proved by integrating Eq. j981 and with the two changes of variable / = — y and t = —u. Using in the last expression the integral 
representation of the Dirac delta, we obtain the desidered result 



f^Hi(y)dy = 7l. 



(100) 



Le us consider Hi (F) = Hi {XY) and study the large F asymptotic behavior of Hi. Supposing small viscosity, Ai<C 1 and X 2 {k 2 + l 2 ) <C 1, and by 
applying Eq. (96} and i97\ to Hi , we obtain for the real part 



while the behavior of the imaginary part is 



XY A 



and 91 \H X {Y)] 



1 



X(k 2 + l 2 ) 

F 2 



(F)] happei 

By performing the change of variables F = y/X, we thus obtain the asymptotic behavior of Hi 



y»i XY ' 

Observe that the crossover in the asymptotic behavior of 9? \H X (F)l happens for F S> 1, so that the asymptotic expansion performed is justified. 



2A3 

" y 4 



and 



for the real part, and 



1 



A 3 (>t 2 + / 2 ) 



(101) 



(102) 



(103) 



(104) 



S[H A (y)] . , 

y»l y 
for the imaginary part. 

A very similar reasoning applied to the function Gi denned in Eq. (67) leads us to the conclusion that 

72 3 

*[<a(y)] ~,-=t- ( 105 ) 

while the leading order behavior of the imaginary part remains the one found for H x . We can thus conclude that H x (y) = Gi{v) + ff{X) for every 
value of v. 



Combination of a viscous and a linear friction regularization. We can now also discuss the regularization of Plemelj formula that cor- 
responds to the combined effect of viscosity and linear friction, by a simple generalization of the argument given above. We start observing that the 
stationary vorticity-vorticity correlations is 



-* 2 ,V1 -yi) =e(k 2 + l 2 )e ,k ^-^e U ^->^ j~ e -^ yi -y 2 )r e -2 v{ kH^r + 2v s k^-^sW g -2a, ^ + c Q 
we have thus to study the following oscillating integral (obtained from the one contained in the above formula by a change of variable) 

Ay JO 



dri, 



(106) 



(107) 



where X v = 



'2v)V3 
ks J 



is the length scale associated with the viscous regularization and X a = ~3 is the length scale associated with the Rayleigh friction. 
The discussion in the following restricts to the real part of Hi v i a because the decay of the imaginary part is the same for all the regularizations. 

Introducing the function Hi v v_ (F) = H Xv i a (X V Y) and performing a reasoning very similar to the one of the previous section, we deduce that the 
real part of H x j. (F) is given by 



X 2 {k 2 + l 2 ) + ^ (x 2 (k 2 + l 2 ) + kA\6Xll{k 2 +l 2 ) + 6X a l + 2 / j 
\ -> h0 i 



Y 2 



F 4 



¥ 6 



(108) 
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We thus see that three different regimes are possible in the limit of small Rayleigh friction (a/s <S 1) and small viscosity (X v l -C 1): (i) the viscosity is 
negligible compared to the Rayleigh friction: X a /X v 2> 1; (ii) the Rayleigh friction is negligible compared to the viscosity X a /X v <C X 2 (k 2 + l 2 ) < 1 
and (Hi) and intermediate case l» X a /X v » X 2 (k 2 + 1 2 ). 

Case (r'): X a /X v 3> 1. We see that equation U08I reduces to 



and thus 



1 

r- 



(109) 



Ay M 



Ml 



and 



(110) 



We thus obtain that, even if the viscosity is negligible with respect to the Rayleigh friction, the behavior 1 jy 2 is modified by the presence of a small 
viscosity in 1 /y 4 for a small enough y, corresponding to X v -C y >C X a . 

Case (ii): X a /X v <C X 2 (k 2 + 1 2 ) <C 1. We see that equation 1108) reduces to the case we have already discussed in the previous section where the 
Rayleigh friction was not present, see Eq. 1 1103) . The only difference with respect to that case is that the crossover point between the 1 /y 4 and 1 jy 2 
behavior is slightly shifted by the presence of a small Rayleigh friction. 

Case (Hi): X 2 (k 2 + I 2 ) <C X a /X v <C 1. With a similar reasoning as before, we obtain 



2X[ 

' v4 



and 



(in) 



We thus see in this case the appearance of the new length scale y corresponding to the crossover between the 1/y decay given by the regularization 
due to the viscosity and the 1 jy 2 decay given by the Rayleigh friction. 



C Convergence of Reynolds stress divergence on stationary points 

We prove in this Appendix that the divergence of the Reynolds stress \(vjm <Bw) (y)J converges even when evaluated at points y = y c , where y c are 

stationary for the zonal flow U: U'(y c ) = 0. Even more genetically, we prove that E{/ KQ)Wi(ri ) f))vm^ (r2,f)))J converges, no matter for which values of 
rj and rj. 

In this case, to the best of our knowledge, no theoretical results are available in literature. However, strong numerical evidences are reported in (7], 
permit to infer relations analogous to equations ( 1701 . J71I . )72t and j73l but for y = y c . The main differences in this case are: (i) the algebraic decay of 
the velocity and of the stream-function depends on the parity of the initial datum; (ii) the exponent of the algebraic decay for y = y c is different. 

For initial conditions even with respect to y — > —y, the numerical results in 1 7 1 can be summarized by replacing Eq. I70K {1 It . 1721 and 1731 by 

1^,01^7 (H2) 

W^A.^i (113) 

| V tW)|,^. (H4) 

Moreover, the asymptotic profile <Bj°(v) appearing in Eq. j70l . )7U . )72l and )73t is equivalent to (y — y c ) 2 in the vicinity of y c . We observe that these 
results are about the asymptotic behavior of the modulus of vorticity and velocity and that the oscillating behavior of these quantities is not known. This 
fact does not allow to obtain estimates for the asymptotic behavior of derivatives of the vorticity and of the velocity. 

For initial conditions odd with respect to y — > — y,the numerical results in |!7] can be summarized by replacing Eq. I70K j71K i72\ and {73} by 

\w k (y c j)\ ~ ^ (115) 



1 

t Jf2 



TCy«0l~_^ (H6) 



i 



,*M^.^- (117) 

Again, the asymptotic profile <5"(y) appearing in Eq. <70t . <7H . 172) and 173) vanishes as (y — y c ) 2 in the vicinity ofy c , 

Using these results, we can prove that the divergence of Reynolds stress at y = y c converges. It is useful here to consider defined, as in section 
E3 by g k ,(xi,X2,yi,y 2 ,t) = e ik ^-^ g k ,(y u y 2 ,t) +C.C. with 

^+<^„+L&,=e"<>'.-«). (118) 
Moreover, because the Orr mechanism works differently for initial conditions with different parity, it is useful to decompose g k l as follow 

#« = gkl.cc-igkl,cs + igkl,sc + gkl.ss (119) 

where gkicc solves the Lyapunov equation 

dgkl.cc .0(1). .0(2) _ . , ... m 

d; +L u k gkl.cc+^u^kSkl.cc = cos(/yi)cos(/y 2 ) ! (120) 

and g k i xs , gki,sc, gkLss solve the same Lyapunov equation with the right hand side replaced, respectively, by cos(/yi ) sin(/\'2), sin(/yi)cos(/v2) and 
sin(/yi ) sin(/y2). The computation needed to show the convergence of the divergence of the Reynolds stress obtained as linear transforms of the four 
functions gu.cc, gkl,cs, gkl,sc and gkl,ss are very similar. We thus report in detail only the one corresponding to gkl t cc- 
With the notation gkl.cc = A? hki iCC and using Eq. j51t . we can write 



hkucc{yuyi)= j *iS*(yi,*i)v^(y 2 ,fi), (121) 
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where <% is the solution to the deterministic evolution d, -\-lS, k and <%(yi .0) = cos(/yi ) and v_ k (y2,0) = cos(/v2)- The convergence of the this integral 
has already been established in the case y\ ^ y c and yi ^ y c m section 15.2.3 1 

For any value of yi and y%, we see from equations l !70t , 11 121 and i ll 14t that the absolute value of the integrand of Eq. d 1 2 1 i decays to zero at 
least as fast as 1 jt 1 . The convergence of hucc i s mus ensured for all the values of yj andy2- 

We leave to the reader the very similar computation needed to show that also the Reynolds stresses obtained from g^ cs , gkl.sc and gu.ss converge. 
The slowest decay of the integrand of expressions like Eq. jl 2 1 1 which is encountered is l/t 3 ' 2 , still sufficient to ensure the convergence of h k \. 

We conclude stressing that the results presented in this appendix rely on conclusions drawn from numerical results reported in 
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